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1. 


SUMMARY 


This  SBIR  Phase  II  Project  has  produced  a  Unified  Flow  Solver  (UFS)  for  Rarefied  and  Continuum 
Gas  Flows  by  collaborative  efforts  of  CFD  Research  Corporation  and  Dorodnizyn  Computing  Center 
of  the  Russian  Academy  of  Sciences.  The  UFS  separates  the  rarefied  and  continuum  flow  domains 
and  selects  appropriate  solvers  to  combine  the  efficiency  of  continuum  models  with  the  accuracy  of 
kinetic  models.  The  Direct  Numerical  Solution  (DNS)  of  the  Boltzmann  Transport  Equation  (BTE)  is 
used  in  rarefied  regions,  while  kinetic  schemes  of  continuum  fluid  dynamics  (CFD)  are  used 
elsewhere.  Using  similar  numerical  methods  for  the  BTE  and  CFD  solvers,  employing  adaprive  mesh 
refinement  technique  and  intelligent  model  selection  algorithms  attain  the  efficiency  and  numerical 
stability  of  the  UFS. 

This  Final  Report  for  Phase  II  work  describes  the  UFS  architecture,  demonstrates  the  capabilities  of 
the  UFS  for  monatomic  gas  flows  and  UFS  extensions  to  molecular  gases  and  reactive  gas  mixtures. 
Using  adaptive  mesh  and  algorithm  refinement  methodology  enables  easy  coupling  of  the  kinetic  and 
continuum  models  within  a  hybrid  code.  Using  kinetic  schemes  for  the  continuum  equations  further 
facilitates  the  coupling.  Several  criteria  for  the  continuum  breakdown  and  domain  decomposition  into 
kinetic  and  continuum  subdomains  have  been  tested.  The  UFS  code  was  parallelized  using  dynamic 
domain  decomposition  and  dynamic  load  balancing  among  multiple  processors  and  tested  at  several 
multi-processor  systems. 

We  have  demonstrated  the  UFS  capabilities  for  several  one-component  gas  flows  and  have  confirmed 
that  hybrid  method  can  result  in  significant  savings  by  limiting  expensive  kinetic  solutions  only  to  the 
regions  where  they  are  needed.  During  the  simulation  process,  the  UFS  can  automatically  introduce  or 
remove  kinetic  patches  to  maximize  accuracy  and  efficiency  of  simulations.  We  have  demonstrated 
the  UFS  capabilities  for  a  wide  range  of  applications  from  hypersonic  external  flows  to  low  speed 
(subsonic)  flows  driven  by  temperature  gradients. 

We  have  extended  the  UFS  to  molecular  gases  with  rotationally  and  vibrationally  degrees  of  freedom 
and  to  multi-component  reactive  gas  mixtures.  Initial  testing  for  mixtures  of  molecular  gases  has 
been  performed.  With  further  development,  the  UFS  can  offer  an  efficient  solution  to  practical 
problems  of  polyatomic  gas  mixtures  of  different  degrees  of  rarefaction. 

The  results  of  this  project  have  been  presented  at  several  conferences  and  published  in  two  journal 
articles.  The  UFS  User  Manual  and  several  Tutorials  have  been  prepared  and  can  be  delivered  to 
potentials  users. 
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2. 


INTRODUCTION 


A  variety  of  gas  flow  problems  are  characterized  by  large  variations  in  the  density  and  other 
macroscopic  characteristics  of  the  gas.  Examples  of  such  problems  include  a  gas  in  the  gravitational 
field  of  a  planet,  where  the  density  may  vary  by  several  orders  of  magnitude  across  the  atmosphere, 
and  a  gas  rotating  at  high  speed  in  a  spinning  cylinder,  with  large  density  variations  in  the  radial 
direction.  Another  class  of  flows  is  characterized  by  the  presence  of  layers  of  relatively  small  extent 
that  are  embedded  in  larger  (near-continuum)  flow  regions.  In  these  layers,  which  are  associated,  for 
example,  with  shock  waves,  contact  discontinuities,  or  shear  or  boundary  layers,  the  state  of  the  flow 
changes  drastically  over  a  relatively  small  distance.  For  space  flight  problems,  the  reduced  density  of 
the  ambient  gas  at  high  altitudes  increases  the  mean  free  path  to  such  an  extent  that  it  becomes 
comparable  to  the  dimensions  of  the  space  vehicle.  During  the  reentry  phase,  such  a  vehicle 
encounters  different  flow  regimes.  At  altitudes  of  90  km  and  above,  the  Knudsen  number  is  large, 
corresponding  to  a  rarefied-gas  regime.  At  altitudes  below  70  km,  the  Knudsen  number  is  smaller  and 
the  flow  is  well  approximated  by  the  continuum  model  using  the  Euler  or  Navier  Stokes  equations. 
For  intermediate  altitudes,  in  the  transitional  regime,  the  continuum  flow  equations  cease  to  be  valid 
in  the  boundary  layers,  but  remain  adequate  in  the  far- field  flow. 

Whenever  the  characteristic  length  of  the  system  (or  the  distance  over  which  a  gas  state  changes 
appreciably)  becomes  comparable  with  the  mean  free  path,  the  classical  continuum  hydrodynamic 
description  of  gas  flow  becomes  invalid.  The  non-equilibrium  kinetic  state  that  exists  in  those 
circumstances  requires  a  full  kinetic  treatment,  by  solving  the  Boltzmann  transport  equation. 
Simulation  of  rarefied  gas  flows  remains  a  challenging  task.  The  rarefied  gas  dynamics  is  the 
synthesis  of  several  fundamental  problems  such  as  molecular  collision  dynamics  and  energy  transfer 
phenomena  in  collisions,  gas-surface  interactions,  condensation  and  evaporation,  plume  and 
expansion  flows,  and  many  others.  All  these  problems  are  related  to  practical  issues  that  can  be 
conventionally  divided  into  two  groups.  The  first  group  covers  the  questions  related  to  hypersonic 
flight  of  vehicles  at  high  altitudes  (mainly  external  flows).  The  second  group  is  mainly  represented  by 
problems  that  involve  material  processing,  micromechanical  devices,  and  microelectronics  (mainly 
internal  flows  with  length  scales  that  are  comparable  with  the  mean  free  path).  Substantial  difficulties 
arising  in  studies  of  both  groups  of  problems  are  caused  not  only  by  the  rarefaction  of  the  gas  but  also 
by  the  presence  of  chemical  reactions. 

Currently,  there  are  several  numerical  approaches  for  solving  problems  of  rarefied  gas  dynamics,  and 
the  choice  of  a  particular  approach  depends  usually  on  the  flow  rarefaction,  the  dimensions  of  the 
problem,  and  the  presence  of  real  gas  effects.  For  weakly  rarefied  flows,  it  is  usually  adequate  to 
account  for  the  effects  of  rarefaction  by  using  slip-velocity  and  temperature-jump  boundary 
conditions  within  the  standard  Navier-Stokes  (NS)  equations.  The  NS  equations  can  be  derived  from 
the  Boltzmann  equation  by  assuming  small  deviations  of  the  distribution  function  from  the 
equilibrium  distribution.  Therefore,  the  NS  equations  become  unsuitable  for  studying  rarefied  flows 
with  finite  Knudsen  numbers  where  the  distribution  function  deviates  appreciably  from  the 
equilibrium  distribution.  The  extension  of  continuum  models  to  rarefied  flows  by  using  the  Burnett 
equations  based  on  the  second  term  of  the  Chapman-Enskog  expansion  has  difficulties  related  to  the 
correct  formulation  of  the  boundary  conditions,  and  the  linear  instability  of  these  equations  to  short¬ 
wave  disturbances. 

For  strongly  non-equilibrium  flows,  one  should  consider  using  the  Boltzmann  equation  and  operating 
within  the  framework  of  the  kinetic  approach.  The  Boltzmann  equation  is  a  nonlinear  integral- 
differential  equation  for  a  one-particle  distribution  function  '.  The  quadratic  nonlinearity  of  the 
collision  integral,  the  dependence  of  the  integrand  function  on  post-collision  velocities,  and  the  high 
multiplicity  (equal  to  five)  of  the  integration  are  the  main  reasons  for  the  mathematical  complexity  of 
the  Boltzmann  equation,  which  distinguish  it  from  other  known  gas  dynamics  equations  and 
complicates  the  use  of  classical  numerical  methods  for  its  solution.  In  the  transitional  regime,  the 
continuum  approach  is  typically  not  adequate  everywhere,  but  kinetic  simulations  are  too  expensive 
because  traditionally  they  require  the  size  of  the  computational  cells  to  be  comparable  with  the  local 
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mean  free  path. 


The  two  main  approaches  that  have  been  used  for  numerical  solution  of  the  Boltzmann  equation  for 
practical  applications  are:  Direct  Numerical  Solution  (DNS)  and  Direct  Simulation  Monte  Carlo 
(DSMC).  Historically,  DNS  is  one  of  the  first  methods  of  solving  the  Boltzmann  equation 2.  It  usually 
consists  of  two  basic  steps:  the  evaluation  of  the  collision  integral  and  the  numerical  integration  of  the 
differential  part  of  the  Boltzmann  equation.  The  major  advantages  of  DNS  are  uniform  accuracy  of 
computing  both  the  low-  and  high-density  regions  of  the  flow,  and  an  easy  and  effective 
parallelization  of  the  computational  code.  The  limitations  of  the  method  are  mainly  related  to 
dependence  of  the  computational  cost  on  the  dimensionality  of  the  problem  and  the  resulting  high 
computational  cost  for  modeling  three-dimensional  problems  3. 

The  DSMC  algorithm  is  based  on  the  assumption  that  a  small  number  of  representative 
“computational  particles”  can  accurately  capture  the  dynamics  of  the  dilute  gas  described  by  the 
Boltzmann  equation  4.  The  DSMC  method  also  splits  the  motion  of  particles  into  two  sequential 
stages:  free-molecular  advection  and  collision  relaxation.  Implementation  of  the  DSMC  method  also 
requires  dividing  the  computational  domain  into  a  grid  of  computational  cells.  The  size  of  these  cells 
should  be  sufficiently  small  so  that  the  change  in  gas  dynamic  properties  across  each  cell  is  small. 
This  size  is  usually  selected  as  the  minimum  value  of  the  mean  collision  time  and  the  mean  residence 
time,  so  that  the  molecules  do  not  cross  more  than  one  cell  during  one  time  step. 

The  DSMC  method  has  become  de  facto  the  main  computational  tool  for  studies  of  complex, 
multidimensional  rarefied  flows  5,6.  This  is  primarily  because  of  its  relative  simplicity,  the  possibility 
of  using  various  models  of  gas  particle  interactions  and  chemical  reactions  without  substantial 
complications  in  the  computational  algorithm,  and  possibility  of  effective  use  of  parallel  computers. 
Compared  to  DNS  approaches,  the  computational  cost  of  DSMC  methods  is  of  the  order  of  the 
number  of  particles.  However,  the  method  becomes  too  expensive  for  near-continuum  flows  because 
its  accuracy  depends  on  the  resolution  of  the  collisional  length  and  the  collisional  time  scale.  Also, 
owing  to  its  statistical  nature,  DSMC  methods  often  yield  low-accuracy  and  noisy  results  relative  to 
DNS  methods,  and  their  convergence  in  general  is  quite  low. 

Aristov  has  presented  a  comparison  of  DNS  and  DSMC  methods  in  a  recent  book 3.  It  was  shown  that 
errors  from  the  use  of  a  uniform  velocity  grid  can  be  smaller  than  statistical  errors  in  DSMC.  The  use 
of  DSMC  can  lead  to  a  small  number  of  particles  representing  “tails”  of  the  distribution  function, 
giving  rise  to  appreciable  fluctuations,  for  example,  in  chemical  reaction  rates.  DNS  in  principle  does 
not  have  this  disadvantage.  Implicit  schemes  in  DNS  can  be  used  to  study  flows  at  small  Knudsen 
numbers  where  the  explicit  schemes  are  rather  slow  in  convergence.  In  the  traditional  DSMC,  explicit 
solution  schemes  are  typically  used.  For  DNS  one  can  use  the  difference  schemes  of  the  second  (or 
higher)  order  of  accuracy,  whereas  accomplishing  higher  order  schemes  with  DSMC  is  not  simple. 
Another  advantage  concerns  the  multiprocessor  computations:  due  to  uniformity  of  grids  in  velocity 
space,  the  parallelization  of  DNS  is  simple  and  effective. 

Efficient  solution  of  the  Boltzmann  equation  can  be  obtained  using  a  spatially  non-uniform  grid  in 
physical  space.  In  the  Knudsen  layers,  one  can  use  a  grid  with  cell  sizes  smaller  than  the  characteristic 
mean  free  path.  In  the  continuum  flow  regions,  one  can  use  the  cell  sizes  much  larger  than  the  mean 
free  path.  In  this  respect,  using  DNS  offers  a  great  advantage  in  comparison  with  the  DSMC  methods, 
which  as  a  rule  have  to  use  computational  cells  smaller  than  the  mean  free  path,  even  in  the 
continuum  flow  regions.  Our  experience  in  using  DNS  for  two-dimensional  and  three-dimensional 
problems  allows  us  to  consider  it  as  a  competitive  method  for  simulating  multi-dimensional  gas  flows 
over  a  wide  range  of  Knudsen  numbers. 

For  several  classes  of  problems,  DNS  is  strongly  preferable  to  DSMC.  The  first  class  of  these 
problems  involves  subsonic  gas  flows.  Using  DSMC  for  simulations  of  subsonic  flows  is  very 
expensive  because  any  perturbations  propagate  with  the  speed  of  sound  (-300  m/s),  whereas  the  gas 
flow  velocity  could  be  much  smaller  (of  the  order  of,  say,  1-10  m/s).  Using  DNS  allows  one  to 
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separate  the  equilibrium  part  and  to  obtain  the  solution  more  efficiently.  The  second  class  of  problems 
involves  gas  mixtures  with  small  fractions  of  some  species.  The  DSMC  method  would  give  large 
fluctuations  in  the  calculation  of  the  macroscopic  properties  of  these  species.  Finally,  DNS  seems 
preferable  for  the  development  of  hybrid  solvers  covering  a  wide  range  of  Kn  numbers  because  both 
kinetic  and  continuum  parts  are  described  in  a  similar  fashion  in  terms  of  partial  differential  equations 
and  finite  volume  numerical  schemes. 

Looking  beyond  computational  methods  for  rarefied  flow  problems,  an  even  greater  challenge  is  the 
development  of  efficient  schemes  for  practical  problems  in  which  the  local  Knudsen  number  varies  by 
several  orders  of  magnitude  within  a  computational  domain,  such  that  regions  of  continuum  and 
rarefied  flows  simultaneously  exist  in  the  domain.  In  such  problems,  the  kinetic  treatment  is  only 
necessary  in  the  rarefied  parts  of  the  domain,  and  a  continuum  approach  is  adequate  for  the  other 
parts.  Since  the  kinetic  treatment  is  much  more  expensive  than  the  continuum  one,  one  should  use  the 
former  only  when  really  needed.  That  is  why  numerical  techniques  based  on  the  combined  use  of 
continuum  and  kinetic  models  are  currently  being  actively  developed. 

In  seeking  a  “unified”  solution  technique,  one  is  confronted  with  several  problems.  First,  one  has  to 
formulate  the  conditions  under  which  a  hydrodynamic  description  (which  is  valid  for  Kn«l)  is 
tolerable.  Second,  one  has  to  select  a  proper  way  to  describe  strongly  non-equilibrium  flows.  In  the 
extreme  case  of  free  molecular  flow,  with  Kn»l,  one  deals  with  a  simpler  problem.  It  is  the 
description  of  the  transition  regime,  with  Kn~l,  and  its  coupling  to  the  continuum  regime,  with 
Kn«l,  which  becomes  the  main  challenge.  Different  methods  of  combining  and  coupling  kinetic  and 
continuum  models  that  have  been  developed  or  proposed  to  date  can  be  classified  into  three 
categories: 

•  Domain  decomposition  in  physical  space.  In  this  category,  the  computational  domain  is 
decomposed  into  kinetic  and  continuum  sub-domains  using  certain  criteria  7,8. 

•  Domain  decomposition  in  velocity  space.  In  this  category,  one  performs  decomposition  in 
velocity  space  to  describe  differently  fast  particles  and  slow  particles  9. 

•  Hybrid  models.  In  this  category,  one  solves  both  the  kinetic  and  the  continuum  equations  in 
the  entire  domain  and  uses  the  distribution  function  to  compute  transport  coefficients  for  the 
fluid  equations  10, 11 . 

The  goal  of  this  Project  is  to  develop  a  Unified  Flow  Solver  (UFS)  for  simulation  of  gas  flows  across 
the  entire  range  of  Kn  numbers  from  the  free-molecular  regime,  to  the  continuum  regime.  In 
particular,  we  develop  a  hybrid  code  that  switches  automatically  from  a  continuum  fluid  dynamic 
solver  to  a  Boltzmann  solver  and  vice  versa.  The  feasibility  of  such  an  approach  has  been 
demonstrated  in  several  recent  papers  12,13,14.  When  coupling  the  Boltzmann  and  continuum  equations 
via  domain  decomposition,  two  problems  have  been  identified:  proper  criteria  for  domain 
decomposition  into  kinetic  and  fluid  regions  and  suitable  numerical  algorithms  for  coupling  the 
different  equations.  Most  works  decompose  the  domain  a-priori.  They  assume  a  Boltzmann  domain  in 
the  vicinity  of  boundaries  and  a  fluid  domain  away  from  the  boundaries.  However,  during  a 
simulation,  the  different  domains  do  not  remain  fixed,  and  cannot  generally  be  accurately  predicted  a- 

priori.  Therefore,  one  needs  criteria  for  automatic  domain  decomposition  such  as  those  proposed  in  15, 

16 

The  uniqueness  of  our  approach  consists  in  using  DNS  for  simulations  of  the  kinetic  domain.  Using 
DNS  instead  of  DSMC  offers  several  advantages  for  building  a  hybrid  code  for  multi-scale 
atomistic/continuum  simulations  of  gas  flows.  One  of  the  major  problems  with  coupling  DSMC  and 
continuum  solvers  is  related  to  strong  fluctuations  of  the  moments  calculated  from  the  DSMC's 
velocity  distribution  functions.  The  problem  of  connecting  the  two  regions  is  generally  overshadowed 
by  rather  severe  stability  problems  when  DSMC  data  are  handed  over  to  the  NS  solver  at  the  interface 
of  the  kinetic  and  continuum  domains.  These  fluctuations  result  in  irregular  boundaries.  Using  DNS 
allows  much  more  manageable  interactions  between  the  continuum  and  kinetic  solvers. 
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The  structure  of  this  Report  is  as  follows.  Section  3  describes  the  UFS  architecture  and  key 
components.  Section  4  is  devoted  to  UFS  validation  and  demonstration  for  monatomic  gas  flows. 
Section  5  describes  UFS  extensions  to  molecular  gases  and  gas  mixtures.  Summary  of  the  work  is 
given  in  Section  6.  The  list  of  publications  and  presentations  about  UFS  is  listed  in  Section  7. 
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3. 


UFS  ARCHITECTURE 


The  key  components  of  the  UFS  are  shown  in  Figure  1.  The  main  component  is  a  Boltzmann  solver. 
For  reasons  described  above,  the  DNS  method  has  been  selected  for  the  solution  of  the  Boltzmann 
equation  implemented  in  this  project.  Another  component  of  the  UFS  is  a  continuum  (Computational 
Fluid  Dynamics)  solver.  It  is  preferable  for  such  a  solver  to  use  numerical  algorithms  similar  to  the 
Boltzmann  solver.  From  this  point  of  view,  the  recently  introduced  kinetic  CFD  schemes  are  very 
attractive  17,18  and  have  been  selected  in  this  work.  The  remaining  components  of  the  UFS  include 
criteria  for  domain  decomposition  into  kinetic  and  continuum  parts  and  coupling  algorithms. 

The  open  source  Gerris  Flow  Solver  (GFS)  19  was  selected  as  a  framework  for  the  UFS.  The  original 
GFS  code  contained  a  binary  tree-based  incompressible  flow  solver  with  a  dynamically  adaptive  grid 
and  support  of  complex  boundaries  20.  The  semi- structured  quadtree/octree  meshes  offer  a  good 
compromise  between  the  flexibility  of  unstructured  meshes  and  the  computational  efficiency 
of  structured  meshes.  Using  the  GFS  framework,  we  have  added  all  the  UFS  components,  first  for  a 
single  component  gas,  then  for  gas  mixtures  and  molecular  gases  with  internal  degrees  of  freedom  of 
the  molecules.  The  UFS  generates  Cartesian  mesh  around  embedded  solid  boundaries  (defined 
through  standard  files),  performs  dynamic  adaptation  of  the  mesh  to  the  solution  and  geometry, 
detects  kinetic  and  continuum  domains  and  selects  appropriate  solvers  based  on  continuum 
breakdown  criteria.  Below  we  describe  the  key  components  of  the  UFS  in  detail. 


Figure  1  Key  UFS  components 


3.1.  Boltzmann  Solver 

The  Boltzmann  transport  equation  (BTE)  describes  the  evolution  of  a  particle  distribution  function /in 
a  six-dimensional  phase  space  21 


/  +  vr.(5D  =  /(/,/)  (1) 

Ot 

Here  r  is  a  position  vector  in  physical  space,  ^  is  the  velocity  vector,  and  t  is  time.  The  right  hand 
side  of  Eq.  (1)  contains  an  integral  operator  describing  binary  collisions  among  particles.  For  elastic 
collisions  in  a  monatomic  gas,  it  has  the  following  form 
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m = f  dco\  (/(§;)/(§’ )  -  f^)m)Scj{g, =-v/+g.  w 

S2  R 3 


Here  v  is  the  collision  frequency,  G  is  the  inverse  collision  integral,  g  =|  ^  —  %  |  is  the  relative 

velocity  of  the  colliding  particles,  co  is  a  vector  on  a  unit  sphere  S 2  in  velocity  space,  and  dco  is  an 
element  of  the  area  of  the  surface  of  this  sphere,  cr(g ,  %)  is  the  differential  collision  cross  section, 
and  x  is  scattering  angle.  The  post-collision  velocities  (^f,  and  the  pre-collision  velocities 
( \ \  i)  satisfy  the  momentum  and  energy  conservation  laws 


(3) 


This  integral  (2)  can  also  be  written  in  the  form 

m  =  J  da | dbj  (f(Of(^)-mi)m)gbd^  (4) 

0  0  R3 


Here  b  is  the  impact  parameter  (defined  as  the  distance  of  the  closest  approach  of  the  trajectories) 
usually  bounded  by  a  certain  value  bm ,  and  8  is  the  azimuth  impact  angle.  The  scattering  angle 
X(g,b)  depends  on  scattering  potential  of  inter-atomic  interactions.  For  the  Hard  Sphere  (HS) 
molecules  of  diameter  d,  the  scattering  is  isotropic  and  b  =  d  sin  9  where  d  —  {n  —  x)l 2 .  For  the 
Variable  Hard  Sphere  (VHS)  models  frequently  used  in  DSMC  simulations,  the  scattering  is  also 
isotropic  and  gcr  =  Ckgl~4/k ,  where  k  is  the  exponent  in  the  intermolecular  potential.  For  other 
commonly  used  scattering  potentials,  these  relationships  can  be  found  in 22. 

We  have  implemented  the  following  scattering  models:  i)  the  HS  model,  ii)  the  inverse  power 
repulsive  potential,  iii)  the  Lennard-Jones  potential,  iv)  the  Coulomb  potential,  and  v)  the  BGK 
model.  For  2D  simulations,  the  BGK  model  was  implemented  in  a  reduced  form  in  which  the  2D 
velocity  space  with  averaging  in  the  z-direction.  Besides  mentioned  potentials  it  is  possible  to  use  in 
the  future  some  modern  potentials  such  as  the  Tang-Toennies  potentials  23 . 

For  the  numerical  solution  of  Eq.  (1),  a  Cartesian  mesh  in  velocity  space  is  introduced  with  a  cell 
size  and  nodes  ^  .  Using  this  velocity  grid,  Eq.  (1)  is  reduced  to  a  system  of  linear  hyperbolic 
system  of  transport  equations  in  physical  space  with  a  nonlinear  source  term 

^  +  Vr-(W  =  I(fpJp)  (5) 

Introducing  computational  grid  in  physical  space,  we  split  the  solution  of  (6)  into  two  stages:  free 
flow  and  relaxation.  For  the  free  flow,  we  used  an  explicit  finite  volume  numerical  scheme 


r*k  rk-\  rk- 1  rk- 1 

JiP  ~  dip  |  p  Ji+l/2,P  ~  Ji-\l2,p  _  q 

At  p  Ax 


(6) 


Here,  the  star  *  denotes  the  intermediate  level,  /^J2  p  is  the  value  of  the  function  on  the  cell  face. 
i  =  (4,4,4)  is  the  3D  spatial  index,  Ax  is  the  3D  spatial  step  (Ax,  Ay,  Az) ,  and  /?  =  (/?,, [i,. ,  fiz )  is  the 
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velocity  index.  The  calculation  of  the  face  values  of  the  distribution  function  we  use  standard 
interpolation  schemes.  For  the  first-order  scheme, 

=  (fi~l  +  f£V  2  -  signer ^  -  /“)/  2 .  (7) 

The  second-order  scheme  has  also  been  implemented  using  three  options:  i)  no  limiter  option,  ii)  the 
so-called  minmod  limiter  (default),  iii)  Van  Leer  limiter.  For  stationary  problems,  the  first-order 
scheme  in  time  with  time  step  selected  from  appropriate  CFL  criterion  is  found  to  be  adequate. 

The  relaxation  stage  has  the  form 


r*k 


fiP  f i/3  _  ,  *k 


At 


*k  r*k  r*k 

Vip  Jip  +  Vy? 


(8) 


Currently,  we  use  an  explicit  scheme  with  automatic  selection  of  the  time  step.  It  is  also  possible  to 
use  implicit  or  explicit-implicit  schemes  for  the  relaxation  stage  to  increase  the  time  step  that  may  be 
important  for  small  Knudsen  numbers.  It  is  also  possible  to  solve  the  system  (6)  without  splitting  into 
collisionless  flow  and  relaxation  and  use  any  other  scheme  for  solving  a  hyperbolic  system  with  a 
nonlinear  source  term. 

We  use  the  Gerris  framework  19  to  generate  dynamically  adaptive  Cartesian  mesh  in  physical  space. 
The  boundary  conditions  specified  at  the  surface  of  solid  objects  imbedded  in  the  computational 
domain  provide  the  distribution  function  of  the  reflected  particles  as  a  sum  of  diffuse  and  specular 
reflections  with  accommodation  coefficient  a 


f(Z>)  =  afM(n,T)  +  (l-a)fr  . 


(9) 


The  specular  reflection  term  is  fr=f  (<%r)  where  is  the  velocity  of  an  incoming  molecule  towards 
the  boundary,  which  after  specular  reflection  transforms  into  \r  —  \  -  2(^  •  w)w  where  w  is  a  unit 
vector  normal  to  the  boundary.  The  diffuse  reflection  term  contains  Maxwellian  distribution 
fM(n,T)  with  a  zero  mean  velocity,  T  is  the  temperature  of  the  boundary,  and  the  density  n  is 

calculated  to  ensure  zero  particle  flux  at  the  boundary  at  a  given  point. 

At  the  boundaries  of  the  computational  domain,  for  most  of  problems,  the  distribution 
function  can  be  assumed  a  Maxwellian  fM(n,u,T)  with  a  mean  velocity  u  for  (^  •  w)  >  0 .  For  the 
parts  of  boundary  with  (^  •  w)  <  0 ,  the  distribution  function  is  found  as  a  result  of  the  solution. 

The  computational  domain  in  velocity  space  is  selected  as  a  box  in  such  a  way  that  the  values 
of  the  distribution  function  outside  of  the  box  are  negligible.  For  2D  (in  physical  space)  problems, 
half  of  the  box  ( <^z  >  0 )  can  be  used. 

3.1.1  Calculation  of  the  collision  integral 

The  main  problem  in  solving  the  Boltzmann  equation  consists  in  evaluating  the  collision  integral 24 . 
The  calculation  of  the  five-fold  integrals  (2)  or  (5)  represents  a  challenge  with  respect  to  efficiency 
and  precision.  We  used  the  discrete  analog  of  the  collision  integral  having  the  following  properties: 

i)  The  integral  should  be  equal  to  zero  for  a  Maxwellian  distribution,  I(fM  ?  fM  )  =  0  . 

ii)  The  distribution  should  remain  positive  for  all  nodes  in  velocity  space  when  the 
relaxation  problem  (8)  is  solved. 

iii)  For  the  collision  invariants,  =  (1,  <^2)  ,  the  conservation  laws  should  be  satisfied 
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(10) 


\y,I(f,f)<%  =  0 

R 3 

Below,  we  review  different  methods  of  calculating  the  integral  and  explain  the  choice  of  the  method 
used  in  this  paper. 

The  NtN  method 

The  first  type  of  methods  can  be  called  Node  to  Node  (NtN).  This  method  has  been  used  by  Goldstein 
25,  Buet 26,  Rogier  and  Schneider 27 ,  Varghese  28  and  Frolova 2.  To  illustrate  the  essence  of  the  method, 
Figure  2  shows  a  collision  sphere  in  velocity  space.  This  sphere  with  center  =  (^.  +§ .)/ 2  and 

radius  \g\/2  is  wrapped  around  pre-  and  post  collision  velocities  due  to  energy  and  momentum 
conservation,  see  Eqs  (4).  The  NtN  method  takes  into  account  only  those  post-collisional  velocities 
that  fall  exactly  into  nodes  of  the  velocity  grid.  Therefore  all  properties  (1-3)  are  satisfied 
automatically.  The  NtN  method  is  conservative  and  deterministic,  and  requires  no  interpolation  of  the 
velocity  distribution  function.  The  drawback  of  the  method  is  that  only  for  defined  scattering  angles 
the  post-collision  velocities  distributed  over  collisional  sphere  fall  exactly  at  the  velocity 
grid  and  the  method  is  applicable  only  for  VHS-like  models  with  isotropic  scattering.  The  NtN 
method  cannot  be  extended  for  more  general  potentials  of  intermolecular  interactions  and  for  non- 
uniform  grid  in  velocity  space  because  only  selected  post-collisional  velocities  are  used.  Also, 
extensions  to  gas  mixtures  with  arbitrary  ratio  of  the  molecular  mass  is  impossible. 

For  the  numerical  calculation  by  the  NtN  method,  the  integral  (2)  can  be  written  in  a  symmetric  form 

28 


1(0 = 2  J dco J  j (8(0  -o  +  so'-o- 8(0  -o- 8(4 - O) fit) AO )g<r(g, (0)dO0 

4  s2  r3r 3 

(11) 

where  is  the  delta- function.  The  advantage  of  this  representation  for  the  numerical  evaluation  is 
that  direct  and  inverse  collisions  are  treated  symmetrically  and  microscopic  reversibility  is  satisfied. 
For  the  numerical  evaluation  of  the  8-fold  integral  (11)  the  integral  is  written  in  the  form 

1  =  4!rZ  llfifj \ +  (!2) 

{  j  s2 

where  A£  is  the  cell  size  in  velocity  space,  A.  =  S(^  -  £)  -S(^t  -  £) ,  A  j  =  -  Q  -  S(^j  -  %)  , 

and  N  is  the  total  number  of  nodes  in  velocity  space. 

To  evaluate  the  remaining  integral  over  the  unit  sphere  ensuring  the  exact  energy  conservation,  the 
NtN  method  takes  into  account  only  those  post-collisional  velocities  that  fall  into  the  collisional 
sphere  (see  Figure  2) 


J  d(°[Ai  +  A;]^(®)  =  Zwy/[A.'  +  A }]  03) 

where  Mfj  is  the  total  number  of  such  points  (for  each  i  and  j)  and  wijt  are  their  weights.  For  the 

VHS-like  models,  the  number  and  position  of  nodes  on  the  collisional  sphere  can  be  determined  a 
priory  and  the  calculation  of  weights  is  simple,  w~l/M.  For  more  general  potentials,  the  angle 
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between  direct  and  inverse  collisions  is  a  function  of  relative  velocity  and  these  calculations  become 
cumbersome. 

The  NtN  method  is  conservative,  micro-reversibility  of  collisions  is  satisfied,  and  no  interpolation  of 
the  distribution  function  is  required.  For  good  accuracy,  a  rather  dense  mesh  has  to  be  used  in  velocity 
space.  If  all  velocity  nodes  of  the  grid  are  used  the  method  is  deterministic.  The  evaluation  of  integral 
in  this  case  requires  0(N 2)  computations.  To  reduce  this  number  of  computations  it  is  possible  to 
select  collisions  using  Monte  Carlo  method 28 . 


Figure  2  Collision  sphere  and  selection  of  the  collisions  for  the  NtN  method 

Tcheremissine’s  method 

To  generalize  the  NtN  method  for  more  complex  cases  it  is  necessary  to  take  into  account  the  inverse 
collisions  that  do  not  fall  exactly  on  the  nodes  of  the  velocity  grid.  Depending  on  how  the  post 
collision  velocities  are  taken  into  account  it  is  possible  to  obtain  different  schemes  of  integral 
calculation.  In  a  series  of  works  (see  29  and  references  therein),  Tcheremissine  has  developed 
conservative  methods  of  calculating  collision  integral  for  arbitrary  interaction  potential,  using 
interpolation  in  velocity  space.  Dividing  contributions  of  post  collision  points  into  two  parts  and 
accounting  them  in  two  of  the  closest  nodes  (see  Figure  3),  it  is  possible  to  satisfy  conservation  laws 
at  each  collision.  This  method  is  briefly  described  below. 

For  arbitrary  potential  on  intermolecular  interactions,  it  is  more  convenient  to  perform  integration 
over  collision  impact  parameters  (5)  instead  of  integration  over  a  unit  sphere  (2).  The  corresponding 
integral  in  the  symmetric  form  is 

m*)  =  ji  ds)  dbj  |  [Sg  -  SO  +  8<g  -  0  +  Sg  -  SI)  +  5<g  -  %)]i&A)bgd^  (14) 

^  0  0  R3  R3 


where  i(^,^l)  =  f(£')f(£l)  —  f(£)f(£ t)  .  The  velocities  before  collision  (^^i)  are  chosen  in 
integer  nodes  N,  Nx  of  the  grid  and  post  collision  velocities  (^  ,  ^  )  may  not  lie  in  integer  nodes.  To 
obtain  the  mass,  impulse,  and  energy  conservation  in  each  collision  and  satisfy  the  condition 
+  ,  the  value  of  +  is  interpolated  to  the  nearest 
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integer  nodes  t;  ,  ^  using  the  following  interpolation: 


^')  +  ^)  =  (1-^)(K^)  +  ^mi))  +  ^(K4)  +  ^i))-  (15) 

On  a  uniform  grid,  it  is  possible  to  perform  this  interpolation  with  one  coefficient  for  five  functions, 
since  conservation  of  mass  and  impulse  in  this  case  is  satisfied  automatically  due  to  symmetric 
position  of  nodes.  The  coefficient  A  is  found  from  the  equation 

& + <£  f = a  - -  xmu  f + (4, )2  ] + mL  f + (#,,  )2  ]  •  (i6) 

The  weight  coefficients  A  and  l-  A  define  the  contribution 

i)  =  )/(£))  t0  the  dosest  integer  nodes  »  see  Fignre  3- 


Figure  3  Selection  of  post-collision  nodes  for  Tcheremissine ’s  method. 

The  interpolation  of  /(<f)/(^i)  can  be  performed  using  any  interpolation  formula  and  does  not 
influence  the  conservation  laws. 

The  NtCN  method 

Finally,  we  describe  a  method  that  can  be  used  for  arbitrary  interaction  potentials  and  non-uniform 
grid  in  velocity  space.  We  start  from  the  same  symmetrical  form  of  integral  as  in  Tcheremissine’ s 
method. 
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% 


Figure  4  Selection  of  nodes  for  NtCN  method. 


The  procedure  of  calculating  collision  integral  consists  of  the  following  steps: 

•  select  values  of  pre  collision  velocities  ^  for  some  impact  parameters  b ,  £ 

•  determine  post  collision  velocities  £  ,  , 

•  find  the  nodes  closest  to  the  nodes  ^  , 

•  make  an  inverse  collision  with  velocities  ^ ^  for  the  same  impact  parameters  b  ,£  , 

•  calculate  velocities  after  this  inverse  collision  £  ,  ^  , 

•  calculate  contributions  to  the  integral  from  the  direct  7y  =  (/^  -  ff^g^b  ,  and  inverse 

Li  =  (  fa.f/i  -  fm  fi  )smi  b  collisions,  where  the  quantities  fjp  ,fj),  are  found  using  a 
logarithmic  interpolation  to  give  zero  integral  for  a  Maxwellian  distribution. 

•  sum  up  contributions  to  the  direct  and  inverse  integrals  —  v(^) /(^)  +  G(^) , 

This  procedure  of  evaluating  collision  integral  uses  closest  nodes  (NtCN)  for  accounting  inverse 
collisions  and  introduces  errors  of  the  order  of  0{h \  f  -  fM  |  in  conservation  of  mass,  momentum 

and  energy.  In  order  to  eliminate  these  errors,  we  introduce  a  correction  to  the  collision  frequency 
using  the  method  of  least  squares 


where  the  coefficients  at  are  defined  from  the  collision  invariants  (11).  Thus,  this  method  of 
calculating  the  Boltzmann  collision  integral  possesses  all  the  properties  (1-3). 

For  evaluation  of  the  eight-dimensional  integrals,  the  Korobov  sequences  30  are  applied.  In  the  general 
case,  Korobov’s  points  in  a  s-dimensional  hypercube  are  defined  as 


(18) 
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where  p  is  a  prime  number,  aPv  are  precalculated  integer  coefficients,  and  the  brace  denotes  the 

reminder  on  dividing  an  integer  by  an  integer.  The  velocity  grid  points  closest  to  the  selected 
Korobov’s  points  are  taken  as  the  velocity  grid  points.  The  accuracy  of  this  procedure  is  estimated  as 

0((ln  Nc)as  /  N '“)  ,  where  the  exponent  a  >  1  depends  on  the  smoothness  of  the  integrated  function. 

For  a  piecewise-constant  function,  a  =  1 .  The  above  error  is  less  than  the  estimated  error  of  the 
Monte  Carlo  method. 

The  typical  number  of  quasi-random  trials  Nc  in  our  simulations  was  equal  to  34000.  We  have  only 

accounted  for  trials  that  fall  inside  a  sphere  (with  center  and  radius  defined  by  the  characteristic 
parameters  of  the  problem)  under  the  condition  that  inverse  collisions  also  fall  into  this  sphere. 
Depending  on  the  value  of  Nc ,  and  the  number  of  cells  in  velocity  space  N,  different  Korobov’s 
sequences  were  selected. 

Note  that  all  three  methods  described  above  make  it  possible  to  solve  the  BTE  without  splitting  into 
the  stages  of  relaxation  and  free  molecular  flow  and  using  any  other  scheme  of  calculating  a 
hyperbolic  system  with  a  source  term. 

The  search  for  the  best  methods  of  calculating  the  Boltzmann  collision  integral  continues.  Many 
attempts  have  been  explored  31,  among  them  are  polar  discretization  of  the  velocity  space  32, 
smoothing  of  the  collision  spheres  33,  and  smoothing  the  collision  integral  34  in  order  to  incorporate 
more  points  of  the  discrete  velocity  grid. 

3.1.2  Validation  of  the  Boltzmann  solver 

In  this  Section,  we  describe  selected  validation  cases  for  the  Boltzmann  solver 

Shock  Wave  Structure 

The  problem  of  shock  wave  structure  is  ideal  for  testing  accuracy  of  the  numerical  BTE  solution,  in 
particular,  the  accuracy  of  the  nonlinear  collision  term.  We  have  performed  simulation  of  the  shock 
wave  structure  for  different  Mach  numbers  and  compared  our  results  with  with  previous  DSMC 
results  and  with  experimental  data.  The  comparison  with  the  benchmark  results  [18]  obtained  by  the 
conservation  splitting  method  and  by  the  Ohwada’s  method  demonstrated  agreement  for  gas  density 
and  temperature  with  accuracy  of  1%  for  the  HS  model  at  Mach  number  M=3. 

Figure  5  compares  our  computations  with  experiments  for  density  and  temperature  profiles  in  the 
shock  wave  in  rare  gases  for  two  different  Mach  numbers.  The  density  profile  in  Argon  for  M=  3.8  is 
compared  with  an  experiment  by  Alsmeyer  35.  The  temperature  profile  in  Helium  for  M=  1.59  is 
compared  to  the  experiment 36.  The  computations  were  performed  for  the  Lennard-Jones  interaction 
potential.  The  agreement  of  experimental  and  calculated  profiles  indicates  the  high  accuracy  of  the 
Boltzmann  solver. 
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Figure  5.  Normalized  gas  density  at  M=3.8  and  temperature  at  M=1.59  for  a  shock  wave  in  argon 
(solid  lines  -  calculations,  symbols  -  experiment) 


Figure  6  shows  longitudinal  and  transversal  temperatures  Tx  and  Ty  in  comparison  with  DSMC  results 
37.  The  velocity  grid  in  the  Boltzman  solver  was  (24,24,12).  One  can  see  that  the  DNS  and  the  DSMC 
results  are  in  close  agreement.  Figure  7  compares  results  of  the  Hard  Sphere  model  and  the  inverse 
power  of  12  model. 
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Figure  6  Comparison  of  the  UFS  results  with  DSMC  results  for  the  shock  wave  problem 

at  Mach  =  5. 
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Figure  7  Comparison  of  Tx  and  Ty  for  two  different  models  of  the  potential: 

HS  -  hard  sphere  model  and  rl2  -  inverse  power  of  12. 

To  evaluate  efficiency  of  the  Boltzmann  solver  for  hypersonic  flows,  we  have  simulated  the  shock 
wave  structure  at  M=  25  using  the  Boltzmann  solver  for  the  hard  sphere  gas.  The  profiles  of 
macroscopic  parameters  are  shown  in  Figure  8  (top  left).  The  slices  of  the  distribution  function  at 
w=h/2  as  a  function  on  v  and  u  velocity  components  in  3  locations:  Xc  +  0.3/1 ,  Xc  -  X  ,  and  Xc  +  X 

are  shown  in  Figure  8.  The  grid  size  in  the  velocity  space  is  /z=1.235  of  thermal  velocity.  The  second 
slice  corresponds  to  the  normalized  density  nearly  0.56,  hence  to  the  ID  graphics  of  the  Muntz  paper. 
The  computations  are  made  at  about  45,000  velocity  nodes  with  CPU  time  of  about  75  hrs. 
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Figure  8  Shockwave  at  M=2  5.  Macro-parameters  and  velocity  distribution  functions  f(v,u,w=h/ 2)  at 

3  positions  inside  the  shockwave. 

Heat  Transfer  Between  Parallel  Plates 

We  have  performed  testing  of  the  Boltzmann  solver  for  a  2  dimensional  problem  of  heat  transfer 
between  two  parallel  plates  with  nonuniform  temperature  along  the  plate  surfaces.  For  the 
collisionless  gas  flow  between  two  parallel  plates  with  temperatures  Tfx)  =  \  at  y=0.5  and 

r2(x)  =  l-0.5sin(2;rx)  at  y=-0.5,  we  obtained  an  analytical  solution  shown  in  Figure  9.  This 
solution  was  used  to  access  limitations  of  the  discrete  velocity  model. 


{K.*M-0.5.O.g  Oty)-<HOL0.5) 


Figure  9  Velocity  distribution  function  at  different  points. 
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In  order  to  test  the  capabilities  of  the  UFS  for  low  speed  problems  we  have  solved  a  series  of 
problems  of  heat  transfer  between  parallel  plates  using  Boltzmann  and  kinetic  NS  solvers.  In  this 
problem,  two  plates  are  heated  to  a  non-uniform  temperature,  which  is  T  —  l-0.5cos(2;rx) ,  along 

the  x-axis,  which  is  the  axis  of  periodicity.  The  results  of  calculations  using  the  Boltzmann  solver  with 
BGK  collision  model  are  presented  in  Figure  10.  One  can  see  that  there  is  also  well-formed  vortex- 
type  distribution  of  the  velocity  flow  field.  Using  the  Boltzmann  solver  for  Kn  from  0.03  up  to  3,  we 
observe  a  maximum  of  the  flow  velocity  of  about  of  0.007  at  Kn=0.l. 


Figure  10  Results  of  simulations  of  the  heat  transfer  problem  using  the  Boltzmann  solver  for  two  Kn 
numbers:  Kn  =  0.03  and  Kn  =  1.  The  total  velocity  profile  is  plotted,  together  with  velocity  vectors. 


Collisionless  Flow  Around  Cylinder 


Consider  rarefied  gas  flow  over  a  cylinder.  At  the  left-hand  side  boundary  (in  front  of  the  cylinder)  we 
assume  that  the  velocity  distribution  function  of  the  gas  is  of  the  form 


(7» 


,  (£-£/J2  +  (£-F0 0)2  +  £\ 

3/2  eXP( - -F- - )■ 


(19) 


where  n^fJ ^ ,  V x  are  the  number  density,  the  velocity  components,  and  the  gas  temperature  in  the 

free  stream.  Introduce  dimensionless  velocity  according  to  U02  =  2 kT0  /  m  .  On  the  surface  of  the 

cylinder,  we  assume  diffuse  reflection  of  the  gas  molecules,  i.e.  the  distribution  of  the  reflected 
molecules  as 
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fr  =  ( nr/nTrf2  exp  (~(^ +^\ 

The  non-penetration  condition  at  the  surface  results  in  the  relation 

j  fr  0,  £  )(£ '  n)d%r  =  -  J  f(t,x,  £  X4  •  «)</£ , 

^■n>0  %-n<  0 

that  defines  the  distribution  function  of  the  reflected  particles,  fr ,  through  the  known  distribution 

function  of  the  incident  particles,  f .  Assuming  the  temperature  of  the  cylinder,  Tr  is  equal  to  the  gas 

temperature  at  the  stagnation  point  on  the  axis,  and  using  the  condition  of  the  constant  enthalpy  along 
the  flow  stream  line,  we  obtain 

Tr=TJ\+(-f^M2), 
where  %  =5/3  and  M  denotes  the  Mach  number. 

Figure  1 1  compares  the  calculated  heat  flux  over  the  cylinder  surface  for  free  molecular  flow  with  an 
analytical  solution.  It  is  seen  that  the  numerical  and  analytical  results  agree  with  high  accuracy. 


Mach  -  2,  Tw  =  2.28,  Free  Molecular  Flow,  ReftneSolid  =  8 


Figure  11  Force  and  heat  flux  for  free  flow  around  cylinder. 

3.2.  Continuum  Flow  Solvers 

Traditional  numerical  schemes  of  Computational  Fluid  Dynamics  (CFD)  are  based  on  discretization 
of  the  continuum  (Euler  or  Navier- Stokes)  equations.  Kinetic  schemes  differ  from  the  traditional  CFD 
schemes  -  they  use  the  BTE  for  building  numerical  CFD  algorithms  (see  38).  Kinetic  schemes  for  the 
Euler  equations  have  been  proposed  in  39,40,  and  independently  in  41,42-  The  main  idea  of  this  approach 
was  suggested  earlier,  in  43 .  Kinetic  schemes  using  moments  of  the  equilibrium  distribution  function 
were  introduced  by  Deshpande  et  al.  44,  and  later  further  developed  and  improved  in  45,  17,  46,  47. 
Generally,  kinetic  schemes  are  preferable  for  hybrid  codes  since  the  BTE  is  used  as  a  foundation  for 
both  algorithms.  We  have  used  kinetic  schemes  for  the  continuum  equations  to  facilitate  coupling  to 
the  Boltzmann  solver.  The  implementation  of  the  kinetic  schemes  for  the  Euler  and  NS  equations  is 
described  below. 

3.2.1.  Kinetic  Euler  Solver 

Our  kinetic  Euler  scheme  follows  the  EFM  Equilibrium  Flux  Method  by  Pullin  39.  The  main  idea  of 
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this  method  is  illustrated  below  for  a  2D  case.  Consider  Euler  equations  in  the  form 


dh  dF  dG  n 

—  +  —  +  —  =  0, 

at  dx  ay 

where 

h  =  {p,pu,pv,pE}, 

F  =  | pu,p  12  +  pu1 ,  pvu,u(pE  +  /?)}, 
G  =  lpv,puv,p/2  +  pv2,v(pE  +  p)  j. 


(20) 


Here  p  is  the  gas  density,  u  and  v  are  the  mean  gas  velocities  along  the  x  and  y  axes, 
correspondingly,  pE  =  3  / 2pT  +  p(u2  +  v2)  is  energy,  T  is  temperature,  and  p  =  pT  is  gas 
pressure  (dimensionless  units  are  used). 


The  discretization  of  Eq.  (20)  using  the  finite  volume  method  gives 


h,n+'  =h,;-—(F" 

,J  ,J  Ax 


-F 

(+1/2 J  1  (-1/2 J 


4y 


(J+l/2 


(21) 


where  hy  is  the  cell  averaged  value  of  h  at  a  time  tn ,  E]"1/2  .  and  G^J+ll2  are  fluxes  on  cell  faces 

along  x  and  y,  correspondingly.  To  obtain  these  fluxes,  we  calculate  the  integrals  over  the  velocity 
distribution  function 


tn+ 1 

Kmj  =7;  j  { V'£xf(xM/2>t>€)d4dt 

txt  t,  R} 


1 ' 


(22) 


Glj+m=-\J^J{yj+m,t,^dt 


where  y/  denotes  the  collision  invariants  defined  above.  The  required  velocity  distribution  at  the  cell 
faces  is  taken  in  the  form 

f  =  H[^]g,  +  (l-H[^)gr,  (23) 

where  gl  and  gr  are  Maxwellian  distributions  at  the  neighboring  cells 


g'  = 


PM  - r  (£,-<U2j)2H£y-<U2j)2 


exp[- 


rjm 
1 (+1/2 


(24) 


and  H[d;\  is  the  step  function 

pU>ol 

[0,^<0j’ 

By  substituting  the  distribution  function  (23)  into  the  expression  for  the  fluxes  (22)  and  performing 
the  integration,  one  obtains 


ma = 


fu,=  J  i.ws’di*  J  e,yg’df 

4x>0  4<0 

G"j+1,2=  f  %yYgldZ+  j  $yvgrd% 

4,>  0  f,<  0 


(25) 
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For  the  first  order  scheme,  the  values  of  macro-parameters  at  faces  /?"_1/2  j  ,  w"+1/2  j  >  U+1/2  j  •>  ^+1/2  j  are 
calculated  for  the  functions  gl  and  gr  using  the  known  values  of  the  macro-parameters  in  cells 
xij(xi+ij)-  For  the  second  order  scheme,  the  calculation  of  these  macro-parameters  is  performed 
using  standard  methods  of  reconstruction  using  the  values  at  x._,  .x{. ,  x.+1  .  (x.  y.x.+1  j ,  x.+2  . )  cells  and 
corresponding  limiters. 

3.2.2.  Validation  of  the  kinetic  Euler  solver 

We  illustrate  the  kinetic  Euler  scheme  for  a  two-dimensional  transient  simulation  of  an  internal  gas 
flow  in  a  channel  with  a  forward-facing  step  at  M=3.  The  first  and  second  order  numerical  schemes 
have  been  used  with  mesh  refinement  based  on  density  gradient.  The  second-order  scheme  employs 
the  minmod  limiter  or  Van  Leer  Limiter.  Figure  12  and  Figure  13  show  simulation  results  with  a 
sensitivity  parameter  for  mesh  refinement  equal  to  0.025.  The  total  number  of  cells  is  5060,  the 
computational  time  is  1  hour,  and  the  memory  usage  is  80  MB  on  an  AMD  2.4  GHz  desktop. 


Figure  12  Computational  mesh  and  gas  density ,  t=4 


Figure  13  Density  contours,  t=4 

Comparison  of  the  results  presented  in  Figure  13  with  published  data  shows  that  the  kinetic  Euler 
solver  gives  results  close  to  those  of  Ref  48. 

3.2.3.  Kinetic  NS  solver 

Minimizing  the  size  of  the  kinetic  domain  where  the  Boltzmann  equation  is  solved  can  increase  the 
efficiency  of  UFS.  Using  the  Navier  Stokes  (NS)  solver  instead  of  the  Euler  solver  can  expand  the 
size  of  the  continuum  domain.  The  idea  of  our  kinetic  NS  solver  is  a  generalization  of  the  scheme 
used  for  the  kinetic  Euler  solver  and  the  Kinetic  Flux  Vector  Splitting  method  by  Chou  and  Baganoff 
45  with  the  distribution  function  at  cell  faces  taken  from  Xu  17.  Details  are  described  below. 

The  development  of  the  kinetic  NS  solver  is  based  on  the  solution  of  the  BGK  equation 

f-  +  vr-(5o  = 

Ct  T 
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where  the  inter-collision  time  r  =  ju/  p  is  expressed  though  gas  viscosity  ju  and  pressure  p.  The 
integration  along  characteristics  gives 

/(r,S,0  =  -  J  g{r\U)e-^dtx  +  e~'lrf0(r  -  %) , 

To 

where  r1  =  r  -  ^{t  -  tl) .  We  use  directional  splitting  method  to  reduce  the  multi-dimensional 
problem  to  a  set  of  one-dimensional  problems.  For  a  one-dimensional  case,  the  distribution  function 
fo~  =  0)  and  the  Maxwellian  distribution  g(x,£x,£y9£z,t)  on  cell  faces  are 

expressed  as  21 : 


/0  =gl[  1  +  +  ^)](1  - tf[x])  +  gr[l  +  arx - r(ar£  +  ^r)]i/[x],  (26) 

g(x,  t )  =  g0  [1  +  (1  -  i/[x])fl ' ?x  +  #  [x]a  rx  +  At] , 
where  functions  al’r ,al>r ,Al’r ,  A  are  polynomial  functions  in  velocity  space 

a  =  ax+  a2%x  +  a£y  +  a£z  +  a5(£*  +  £y  +  4;) 


with  coefficients  al,r,al,r  expressed  through  gradients  of  the  macro  parameters  in  physical  space. 
Coefficients  AI,r  are  calculated  to  satisfy  conservation  laws: 

\{a,'Zx  +  A,')yad4  =  0, 

where  y/a  are  the  collision  invariants.  The  parameters  of  the  Maxwellian  distribution  g0  are 
calculated  from  the  relation 


J 'goVad%=\  S'Vad^+\  S>ad Z’ 

4>o  4<o 

For  calculation  of  function^,  the  condition  suggested  in  49  is  used 

Jt  J[g(0,£0  -  u=  0 . 

Having  obtained  the  velocity  distribution  function  on  cell  faces,  the  particle  fluxes  on  faces  are 
calculated  by  integration  of  the  velocity  distribution  function  with  the  collision  invariants 

Fa,x+l/2  =  If*. 

Ga,y+ 1/2  =  If*.  f(yt+ 

This  scheme  incorporates  the  non-equilibrium  character  of  the  distribution  function  by  additional 
parameter  TgI,r(al,r <4  +  AI,r )  and  approximates  the  NS  equations  if  Tgl,r(al,r <4  +  AI,r )  «  1 
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3.2.4  Validation  of  the  kinetic  NS  solver 


Figure  14  shows  results  of  simulations  for  gas  flow  around  a  prism  for  M=  5,  angle  of  attack  of  6 
degrees,  using  kinetic  NS  solver.  The  boundary  condition  at  the  body  surface  is  diffusive  reflection 
with  the  temperature  of  the  wall  T=9.  Mesh  adaptation  performed  using  parameter 
8  =  log(ft)  +  log (u)  +  log(v)  where  n  is  the  gas  density,  u  is  longitudinal  velocity  and  v  is 
transversal  velocity.  This  value  of  the  parameter  gives  correct  concentration  of  the  grid  not  only  in  the 
area  of  the  shock  waves  but  also  behind  the  body  that  is  important  for  study  vortex  flows. 


Figure  14  Gas  flow  around  a  prism  at  M=5,  Kn=0.0001.  The  density  (upper  left)  and  velocity 
distribution  (lower  left),  the  computational  grid  (upper  right)  and  velocity  vectors  behind  the  body 

(lower  right). 


The  results  of  calculations  of  low  speed  gas  flow  between  nonuniformly  heated  plates  are  presented  in 
Figure  15.  One  can  see  that  there  is  a  well-formed  vortex  in  the  middle  of  the  simulation  domain.  This 
vortex  is  temperature-driven.  The  flow  velocity  at  larger  Kn  =0.03  is  about  twice  as  large  as  for 
Kn=0.001  and  its  value  of  0.0012,  which  is  in  very  good  agreement  with  the  results  obtained  by  Sone 
38  using  the  linearized  Boltzmann  calculations. 


Figure  15  Results  of  simulations  of  the  heat  transfer  problem  using  the  kinetic  NS  solver  for  two  Kn 
numbers:  Kn  =  0. 001  and  Kn  =  0. 03.  The  total  velocity  profile  is  plotted,  together  with  velocity 

vectors. 
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3.2.5  Prandtl  Correction 


It  is  well  known  that  the  BGK  model  results  in  incorrect  Prandtl  number,  Pr=l.  To  introduce  Pr 
correction,  we  calculated  the  heat  flux  FH  on  the  cell  faces  using  polynomial  interpolation  of  the 
velocity  distribution  function  defined  at  cell  centers.  The  resulting  heat  flux  has  the  form 

fe=fe+(P-\)fh. 

Pr 


This  algorithms  of  Pr  correction  was  tested  for  the  shock  wave  structure  for  different  M  numbers 
(1.5,3,5,10)  and  for  different  temperature  dependence  of  the  viscosity  coefficient  (see  Figure  16).  The 
results  of  calculations  were  compared  with  the  benchmark  calculations  by  Xu  using  the  classical  NS 
solver. 


Shock  Structure  Pr=2/3  Kn=:1  Ma=1Q 


Shock  Structure  Fr=2/3  Kn=1  Ma=5 


Figure  16  Comparison  of  the  kinetic  NS  solver  with  Pr  correction  with  classical  NS  solver  for  M=10 

(left)  and  M=5  (right) 


The  BGK  model  implemented  in  the  new  form  allowed  us  to  implement  the  Shakhov  correction  in  an 
simple  and  efficient  manner.  The  new  BGK  model  with  Shakhov  correction  has  been  benchmarked 
for  the  shockwave  problem  as  different  Mach  numbers.  Figure  17  shows  results  of  comparison  for  the 
heat  flux  between  3  models:  the  BGK  model  with  no  Pr  correction  (Pr  =1),  the  BGK  model  with  Pr 
correction  (Pr  =  2/3)  and  the  full  Boltzmann  calculation.  One  can  see  that  the  BGK  model  with 
Shakhov  correction  reproduces  very  well  the  results  of  full  calculations  and  that  the  results  without  Pr 
correction  differ  significantly  from  the  full  calculations  (see  also  Figure  18). 


Figure  1 7.  Comparison  of  the  heat  flux  profiles  calculated  for  ID  SW  at  Mach  =  2  using  the  full 
Boltzmann  collisional  integral  (HS  model),  the  Shakhov  Model  with  Pr  =  2/3  and  the  BGK  model  with 

Pr  =  1. 
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Figure  18.  Normalized  difference  in  heat  flux  calculated  with  and  without  Pr  correction  for  the  SW 

problem  at  Mach  =  2. 


Figure  19  shows  the  results  of  UFS  calculations  of  SW  structure  at  M=2  with  the  BGK  model  without 
Pr  correction  (Pr  =  1)  and  with  Pr  correction  (Pr  =  2/3). 


0  5  10  15  20  25  30  35  40 

x/? L 


0  5  10  15  20  25  30  35  40 


x/k 


Figure  19.  UFS  calculation  of  SW  structure  for  two  Pr  numbers:  Pr  =  1  (BGK  model)  and  Pr  =  2/3 

(Shakhov  Model). 

The  BGK  model  with  Pr  correction  has  been  further  expanded  to  the  case  of  3T-BGK  model 
(described  below)  describing  molecular  gases  with  internal  degrees  of  freedom  of  the  molecules.  The 
results  of  calculations  of  SW  at  Mach  =  5,  Zr  =  3,  Zv  =  100  and  Pr  =  2/3  are  shown  in  Figure  20. 
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SW:  3T-BGK,  N2,  Mach  =  5,  Zr  =  3,  Zv  =  100,  Pr  =  2/3 


Figure  20.  Calculation  of  ID  SW  using  3T-BGK  Model  with  Shakhov  correction  with  Pr  =  2/3. 

3.3.  Domain  Decomposition  Criteria 

The  main  problem  of  unified  methods  is  to  separate  kinetic  and  continuum  regions.  In  our  solver  the 
adequate  switching  criterion  is  important  because  wrong  domain  decomposition  could  lead  to  non¬ 
positive  distribution  function  when  kinetic  NS  solution  is  coupled  with  the  Boltzmann  solution.  We 
have  used  the  following  switching  criteria: 

Sp=Jp2„  +  p2xy+pl/p.  (27) 


J  Kn  local 


1  dp 

p  dx  ’ 


(28) 


Sgkp=Kn 


'Yr\ 

V  p  J 


1 

+  — 
T 


r a.,\2  r^.\2 


du 

\dx  j 


+ 


dv 

\dyj 


+ 


v5zy 


(29) 


s ns  —  r.n 


'vpY 

V  P  V 


+ 


r  du 
Kdx 


+ 


f  dv^2  ' 


j 


+ 


dw 
\  dz  y 


l 


2  2  2 

U  +  V  +  w 


)  (30) 


where  pxx ,  pyy ,  pzz  are  appropriate  components  of  the  non- equilibrium  stress  tensor,/?  is  the  pressure, 
p  is  density,  T  is  temperature,  u,  v,  w  are  appropriate  component  of  velocity  (all  values  are  given  in 
dimensionless  form),  Kn  is  the  Knudsen  number  of  the  problem  under  consideration  (e.g.,  for  a  flow 
around  a  cylinder  Kn  =  A  /  R  where  A  is  the  mean  free  path  and  R  is  the  radius  of  a  cylinder).  If  S  is 
greater  then  a  threshold  value,  then  the  kinetic  solver  must  be  used.  The  applicability  of  different 
criteria  and  the  ways  to  choose  the  threshold  value  is  currently  being  studied.  It  was  found  that  the 
criterion  SKn  local  gives  correctly  the  non-equilibrium  domain  near  shock  wave  and  behind  the  body 
at  moderate  Knudsen  numbers,  but  at  small  Knudsen  numbers  (Kw«0.1)  non- equilibrium  domain 
behind  the  body  appears  to  be  too  small.  We  found  experimentally  that  criterion  SNS  gives  correct 
kinetic  regions  and  allows  one  to  successfully  couple  the  NS  and  Boltzmann  solvers. 
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We  have  studied  the  influence  of  the  breakdown  parameter  on  the  flow  characteristics  calculated  by 
the  UFS.  Figure  21  shows  an  example  of  the  normal  force  on  the  cylinder  surface  calculated  at  5 
breakdown  parameters  s  =  0.025,  0.05,  0.1,  0.2,  and  0.4  for  supersonic  gas  flow  around  a  cylinder  at 
Mach=3.  One  can  see  that  all  curves  converge  at  small  s  numbers  when  the  Boltzmann  region  grows. 
At  the  same  time,  by  decreasing  the  s  number,  the  computation  time  increases.  Therefore,  for  quick 
results  one  can  use  larger  s  numbers  if  precision  of  the  order  of  10%  is  satisfactory. 


Fnormal 


Figure  21  Normal  force  on  the  cylinder  body  as  a  function  of  angle  for  the  problem  of flow  around 

cylinder  for  Mach= 3. 

3.4.  Coupling  kinetic  and  continuum  solvers 

The  coupling  of  kinetic  Boltzmann  and  continuum  Euler  solvers  consists  of  the  following.  The  Euler 
equations  are  solved  in  the  entire  computational  domain.  The  boundary  conditions  for  the  Boltzmann 
equation  at  the  kinetic/continuum  interface  are  obtained  assuming  Maxwellian  velocity  distribution 
function  in  the  continuum  cells.  In  the  kinetic  domain,  the  moments  are  obtained  from  the  Boltzmann 
solution. 

The  coupling  of  the  Boltzmann  solver  and  the  NS  solver  consists  of  the  following.  On  each  time  step, 
a  continuum  cell  is  considered,  which  is  a  neighbor  to  a  kinetic  (Boltzmann)  cell.  In  this  continuum 
cell,  a  velocity  grid  is  introduced  which  is  identical  to  that  in  the  kinetic  cell.  On  this  velocity  grid,  the 
following  distributions  functions  are  constructed  f0  =  g-z[l  —  T(al^n  +  Alf\  on  each  face  where  is 

the  normal  velocity  to  a  cell  face.  The  parameters  of  the  Maxwellian  distribution  function  gzare 
calculated  using  the  macroparameters  in  the  continuum  cell  and  the  coefficients  of  the  polynomial 
a1  are  calculated  using  the  gradients  of  the  macroparameters  in  the  continuum  and  the  neighboring 
kinetic  cells.  The  coefficients  of  the  polynomial  A1  are  then  calculated  using  the  relationship  of 
conservation  J  gl(al^n  +  Al)y/ad ^  =  0  of  the  moments  on  the  discreet  velocity  grid.  Coupling  NS 

and  Boltzmann  solvers  requires  knowledge  of  the  face  values  of  the  distribution  function,  whereas 
only  cell  values  are  used  for  coupling  with  the  kinetic  Euler  solver.  The  coupling  with  the  kinetic  NS 
requires  that  the  face  values  are  stored  and  transferred  to  the  Boltzmann  solver. 

The  results  of  coupled  NS/Boltzmann  solution  are  shown  in  Figure  22  for  the  ID  shock  wave.  The 
Boltzmann  solver  was  run  with  the  HS  model  for  the  collisional  integral.  One  can  see  that  pure 
Boltzmann  results  are  very  close  to  those  obtained  using  the  coupled  NS-Boltzmann  solution.  Only 
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about  a  third  of  the  simulation  domain  is  simulated  by  the  Boltzmann  solver,  that  gives  significant 
computational  speedup. 


SW:  Boltzmann  and  Boltzmann/NS,  Mach  =  3 


distance  /  lambda 

Figure  22  Results  of  comparison  between  the  pure  Boltzmann  (symbols)  and  NS-Boltzmann 
computations  for  the  shock  wave  at  Mach  =  3.  Also  shown  is  the  kinetic  flag  indicating  the  region 

(flag  =  1)  where  the  Boltzmann  solver  is  run. 

3.5.  Axi-Symmetric  version  of  UFS 

The  axi-symmetric  version  of  the  Boltzmann  solver  (with  BGK  collision  term)  has  been  implemented. 
The  solver  uses  cylindrical  coordinates  (x,r,(p)  in  physical  space  and  cylindrical  coordinates 

in  velocity  space.  The  cylindrical  velocity  components  are  related  to  the 

Cartesian  velocity  components  as 


4  = 

=  4y  cos  (P  +  %z  sin  cp,  (31) 

^  =  -4vsin^+£cos^ 


In  these  coordinates,  the  Boltzmann  equation  (1)  has  the  form 

3,0/)  +  %xdx(rf)  +  ldr(rf )  +  3&  (£/)  -  0^  (££/)  =  rl .  (32) 

where  it  is  assumed  that - /  =  0  .  As  pointed  out  in  50,  for  the  numerical  solution  of  Eq.  (50)  it  is 

dcp 

more  convenient  to  use  cyclic  coordinates  (R,a>)  defined  as:  %r=R  cos  co,fp  -  R  sin  co , 

+  ^  .In  these  coordinates,  Eq.  (32)  has  the  form: 

3,0;/')  +  <Hxdx(rf )  +  R  cos  adr{rf)  -  dj&m&f)  =  rl .  (33) 

We  used  collision  integral  in  the  BGK  form,  /  =  (fm  —f)/z,  where 
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fm  =  n(l/ 7rT)3  2  exp|-[(^  -  ux )2  +  (7? cos co-ur)2  +(7?sin&>)2]/rj  is  the  Maxwellian 

distribution.  The  local  parameters  of  the  Maxwellian  distribution  fm  are  determined  at  each  time  step 
according  to  relations: 

( n ,  nux ,  nur  ,  2?)*  =  J  (1,  gx ,  7?  cos  (<^2  +  7?2  )/  fRd  %xdRd  co , 

£  =  «(|r  +  K2+«2)). 


For  the  numerical  solution,  we  introduce  uniform  mesh  in  velocity  space  {£i9Rj,G)k}  with  steps 

{ A<^,A/?,A&>}.  The  approximation  of  the  convective  terms  dx(rf  ),R  cos  cod r(rf)  is  done  by  the 

standard  procedures  [50].  The  differential  approximation  of  the  term  d a (sin  cof)  must  satisfy 

additional  conditions  to  give  correct  values  to  the  discrete  analogues  of  integrals  from  trigonometric 
functions,  which  is  necessary  to  satisfy  conservation  laws  and  ensure  positive  value  of  the  velocity 
distribution  function  [50].  We  used  the  following  approximation: 


djsmcof) 


1 

2sin(A<y/2) 


[(sin^+i/2)+A+i  +  (sin<qt+i/2)"/t -(smo)k_l/2)+ fk  -  (sin cok_v 2)~  fk_{\ 


where 


—  (a±|a|)  and 


c°k±i/2  =  —  l/2Aty.  As  shown  in  [50],  such  an  approximation 


ensures  that: 

1 .  The  positivity  of /is  preserved. 

2.  The  conservation  laws  of  density  n  ,  impulse  nux  and  energy  E  are  satisfied. 

3.  The  entropy  is  locally  dissipated. 

4.  The  uniform  flows  are  preserved. 


To  derive  kinetic  scheme  for  the  gas  dynamic  equations  in  the  cylindrical  system,  the  kinetic  equation 
(50)  is  integrated  over  velocity  space  with  invariants  y/1  =  (l,^x,^r,(^x  +  <^2  +  <^2))  to  obtain 


d,(rX)  +  %xdx(rY)  +  £dr(rZ)  +  F  =  0 


(34) 


where 

x  =  \fy'd4xd4M9> 

Z  =  \U¥d^d^, 


F  =  J  (0&  (£/)  -  ^  (££  /))  • 

For  a  Maxwellian  distribution  f, 
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X  -  {n,nux,nur,nE} , 

Y  =  ^nux,p / 2  +  nux  ,nuxur,ux(nE  +  /?)j, 
Z  =  {nur,nuxiir,  p  /  2  +  nur2  ,ur(nE  +  p)  | 


F  =  {0,0,/?/2,0}. 

The  numerical  scheme  for  solving  Eq.  (34)  has  the  form 


X  -  X  n 
At 


f  yn  _yn  r  7n  -r  7”  ^ 

1  i+l/2J  1  i-l/2,j  [j+ 1/2^  ij+ 1/2  0-1/2^  ij -1/2  p 


Ax 


A  r 


hj 


y 


(35) 


To  calculate  fluxes  on  cell  faces,  we  used  /  = 


\fUL>  o)] 
l/;,«,<o)J 


it,-*':)2 +&-<' r+i 


l,r\2 


fLr  =  . 
J  M 


-exp[- 


where 


-] ,  and  the  upper  indexes  (l,r)  correspond  to 


(. nTl'r)vl  Tl’r 

the  left  and  right  cell  center  values,  and  is  the  normal  component  of  the  velocity  at  the  face.  This 
way  we  obtained  kinetic  scheme  for  the  Euler  equations  in  the  cylindrical  coordinate  system. 

Together  with  Eq.  (34)  we  used  the  equations  obtained  from  the  non-conservative  form  of  the  kinetic 
equation  (33): 

SX  dYdZF 

- +  —  +  —  +  — =  0,  (36) 

dt  dx  dr  r 


where  F  =  j nur,nuxur,nu2,(E  +  p)ur}  .  As  pointed  out  by  many  researchers,  these  equations  are 
preferable  for  the  description  of  the  solution  near  the  axis  of  symmetry. 

The  axisymmetric  version  of  the  UFS  code  has  been  tested  against  full  3D  calculations  for  a  quarter  of 
a  supersonic  nozzle.  The  results  of  the  tests  are  shown  in  Figure  23. 


distance.,  cm 


)  Pure  N*  /  M  *  I  i  -  3  7xt  0-2  W m3 1  Tr(H  -300  K 


Figure  23.  Comparison  between  2D  axi-symmetric  and  full  3D  calculations.  Shown  are  the  density 
and  translational  and  vibrational  temperature  profiles  along  the  central  line. 
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3.6.  Parallelization 


The  initial  version  of  the  GFS  code  had  possibility  of  parallel  domain  decomposition  using  cubical 
boxes  or  squares  (in  2D).  For  parallel  execution,  the  computational  domain  could  be  subdivided  into 
several  cubical  (square)  sub-domains,  and  a  selected  set  of  these  sub-domains  assigned  to  different 
processors.  Such  decomposition  is  static  (since  it  does  not  change  during  the  computation)  and  rather 
inefficient.  First,  for  the  domain  decomposition  into  the  cubic  sub-domains  is  not  always  possible  to 
achieve  good  load  balance  between  the  processors,  even  for  the  cases  when  the  load  is  known  for  each 
computational  cell,  and  the  computational  grid  is  also  known.  When  adaptive  grid  is  used,  static 
domain  decomposition  can  result  in  large  load  disbalance.  Moreover,  for  the  UFS,  the  computational 
load  in  each  cell  could  wary  by  orders  of  magnitude  depending  on  weather  the  cell  is  kinetic  or 
continuum. 

For  these  reasons,  it  was  decided  to  abandon  the  parallel  option  built  into  original  GFS  and  develop 
new  parallel  capabilities  using  subdomains  of  arbitrary  shape  with  dynamic  load  balancing  between 
the  processors  depending  on  local  grid  refinement  and  different  weight  of  kinetic  and  continuum  cells. 
This  problem  was  subdivided  into  several  stages. 

First,  the  capability  of  performing  computations  in  a  selected  part  of  the  computational  domain  of 
arbitrary  shape  was  implemented  in  UFS.  To  illustrate  details  of  the  implementation,  consider  the 
procedure  of  accessing  different  cells  in  the  GFS  code  illustrated  in  Figure  24.  The  computational  grid 
in  the  GFS  is  generated  by  subsequent  division  of  square  boxes  into  smaller  boxes  with  linear 
dimensions  equal  to  half  of  the  initial  dimension  (left  part  of  Figure  24).  The  procedure  of  creating  the 
computational  mesh  can  be  represented  by  a  tree  (the  right  part  of  Figure  24).  The  root  of  the  tree  (0th 
level)  corresponds  to  initial  cube;  the  first  level  corresponds  to  4  (in  2D)  or  16  (in  3D)  cubes  obtained 
by  division  of  the  initial  cube,  etc.  The  computational  cells  correspond  to  leaves  of  the  tree.  The  order 
of  cell  traversing  is  shown  in  Figure  24  by  dashed  lines  with  arrows.  All  computation  procedures  are 
called  only  for  the  leaves  of  the  tree. 


Figure  24  The  cell  traversing  procedure 

To  perform  computations  only  for  a  part  of  the  domain  (for  instance,  the  sub-domain  shown  by  blue 
color  in  Figure  25),  one  introduces  a  flag  for  each  leaf  to  identify  whether  or  not  the  cell  belongs  to 
the  selected  sub-domain.  If  a  cell  is  a  parent  cell  for  a  set  of  cells,  it  is  flagged  only  if  at  least  one  child 
cell  belongs  to  the  selected  sub-domain.  After  introducing  flags,  the  procedure  of  cell  traversing  is 
modified  in  such  a  way  as  to  visit  only  the  cells  belonging  to  the  selected  sub-domain.  As  shown  on 
the  right  part  of  Figure  25,  only  the  branches  connected  by  solid  lines  are  traversed.  Moreover,  it  is 
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necessary  to  specify  boundary  conditions  at  the  boundary  of  the  selected  sub-domain.  For  this 
purpose,  additional  cells  marked  by  dashed  lines  in  Figure  25  are  used.  The  branches  corresponding  to 
the  boundary  cells  are  shown  by  dashed  lines  on  the  graph.  The  boundary  cells  are  marked  using  a 
different  flag,  and  can  be  traversed  separately  by  the  code. 


■ 

1 

m 

iH 

Figure  25  Simulations  of  a  part  of  the  domain 

The  implementation  of  the  approach  described  above  enabled  us  to  perform  simulations  only  in 
selected  parts  of  the  computational  domain.  An  additional  outcome  from  this  part  of  the  work  is  the 
possibility  to  simulate  one-dimensional  problems  (in  the  original  GFS,  these  problems  had  to  be 
solved  as  two-dimensional  problems). 

The  decision  on  the  necessity  to  perform  load  balancing  is  made  taking  into  account  the  following 
algorithm.  We  consider  the  maximum  and  average  (over  processors)  computational  time  for  each  time 
step  Tmax  and  Ta.  We  introduced  2  coefficients  Afast  and  Asiow.  The  load  must  be  balanced  if 
Tmax>  Taf  1  +Asiow+Afast/Tmax).  The  value  Asiow+Afast/Tmax  corresponds  to  deviation  of  the  parallel 
efficiency  of  the  computational  part  of  the  program  from  unity.  The  first  term  in  this  sum  is  the  main 
term  by  which  we  specify  the  limit  for  the  above  deviation.  The  second  term  is  necessary  for  the  fast 
computations,  where  the  repartitioning  time  is  comparable  with  the  computational  time  and  where  too 
frequent  repartitioning  can  give  rise  to  significant  deceleration  of  the  algorithm.  The  coefficients  A/ast 
and  Ashw  can  be  specified  by  user  in  sim-file,  and  should  be  chosen  different  for  different  computers, 
depending  on  the  computational  performance  and  data  exchange  speed.  This  criterion  for  load 
balancing  proved  effective  for  Euler  and  Navier-Stokes  solvers,  for  Boltzmann  solver,  and  for  unified 
computations  engaging  both  kinetic  and  gas-dynamics  solvers. 

The  next  stage  of  the  UFS  parallelization  consisted  in  static  parallelization  when  different  sub- 
domains  were  assigned  to  different  processors  and  remained  unchanged  during  the  computations. 
After  completion  of  this  stage,  the  possibility  of  dynamic  domain  decomposition  was  implemented. 
The  dynamic  load  balance  between  different  processors  was  achieved  separately  for  kinetic  and 
continuum  equations.  The  procedure  of  domain  decomposition  was  performed  using  space-filling 
curves  (SFC)  as  illustrated  in  Figure  26.  During  sequential  traversing  of  the  cells  by  natural  order,  the 
physical  space  is  filled  with  curves  in  N-order  (Morton  ordering).  After  this  ordering  of  cells,  all  cells 
can  be  considered  as  a  one-dimensional  array.  A  weight  is  assigned  to  each  cell,  which  is  proportional 
to  the  CPU  time  required  to  perform  computations  in  this  cell.  Furthermore,  the  array  modified  with 
corresponding  weights,  is  subdivided  into  sub-arrays  equal  to  the  number  of  processors,  in  such  a  way 
that  the  weight  of  the  sub-arrays  are  approximately  the  same.  This  method  allows  rather  efficient 
domain  decomposition  between  different  processors.  Figure  27  shows  an  example  of  domain 
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decomposition  between  several  processors  for  supersonic  gas  flow  around  a  cylinder.  Figure  28  shows 
an  example  of  dynamic  load  balancing  for  a  3D  problem  of  gas  flow  around  a  Space  Shuttle  Orbiter. 


Figure  26  Space  Filling  Curves 


Figure  27  An  example  of  DLB  between  8  processors  for  the  problem  of flow  around  cylinder.  Each 

processor  ID  is  shown  by  a  different  color. 


Figure  28  An  example  of  DLB  between  8  processors  for  the  problem  of flow  around  the  Space  Shuttle 
Orbiter.  Each  processor  ID  is  shown  by  a  different  color. 
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4. 


UFS  VALIDATION  AND  DEMONSTRATION  FOR  MONATOMIC  GAS  FLOWS 


4.1.  High  Speed  Flows  Around  Blunt  Bodies 

4.1.1.  Flow  around  a  cylinder 

The  flow  around  a  cylinder  at  M=  3  was  studied  in  a  wide  range  of  Kn  numbers  (0.001-10)  using 
different  criteria  of  domain  decomposition  and  the  temperature  of  the  cylinder  surface.  Figure  29 
shows  the  results  of  computations  with  Boltzmann  solver  and  Euler  solver  using  diffuse  reflection 
with  the  wall  temperature  of  T=  4.  The  switching  criterion  is  SNS,  the  breakdown  parameter  is 

Ss  =  0.3  .  The  parameter  of  the  mesh  refinement  is  8  -  log(/?)  +  log (u)  . 


Figure  29  Gas  flow  around  a  cylinder  for  M=3  for  different  Kn  numbers  ( Kn  =  A  /  R  =  5,  1.5,  0.5, 
0.05,  0.005).  On  the  left  side  are  the  density  profiles,  on  the  right  side  are  the  computational  grid  with 

kinetic  (red)  and  continuum  (white)  domains. 

Important  quantities  of  the  gas  flow  interaction  with  solid  bodies  are  the  vector  of  forces  ( Fx, ,  Fy,  Fz ), 
the  normal  force  Fn  (pressure)  and  the  shear  force  or  shear  stress  Ft  and  the  heat  flux  H.  Both  integral 
and  local  quantities  can  be  calculated  by  the  UFS  code.  We  have  tested  this  implementation  for  the 
problem  of  gas  flow  around  a  cylinder  at  M=  3.  The  results  are  shown  in  Figure  30,  Figure  31,  and 
Figure  32  for  different  levels  of  mesh  refinement  near  the  cylinder  surface.  One  can  see  that  the  shear 
force  is  a  relatively  smooth  function  around  the  cylinder.  Although  not  directly  compared  here,  the 
shear  force  appears  to  be  smoother  than  what  is  expected  from  a  cut-cell,  conventional  Navier-Stokes 
formulation  51,  52.  As  noted  in  51, 52, 53,  the  non-smoothness  and  non-orthogonality  at  mesh  refinement 
and  cut-cell  boundaries  of  hierarchically-based,  adaptively  refined  grids  can  introduce  non-positive 
discrete  representations  when  applied  to  solving  the  Navier-Stokes  equations.  This  is  due  to 
representation  of  the  viscous  fluxes  by  higher  order  derivatives  of  the  cell-centered  data  in  a  non¬ 
positive  fashion,  which  for  low  cell  Reynolds  and  Peclet  numbers,  can  cause  at  best,  non-smooth 
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solutions,  and  at  worst,  instabilities  that  may  not  be  damped  by  the  temporal  scheme  51 .  It  seems  the 
approach  implemented  in  UFS  might  overcome  these  deficiencies,  although  more  investigations  are  in 
order. 


Figure  30  Shear  and  normal  force  for  the  problem  of  gas  flow  around  a  cylinder. 


angle,  degree 


Figure  31  Heat  flux  and  Fx  for  the  problem  of  gas  flow  around  a  cylinder 


We  have  further  verified  that  increasing  the  mesh  resolution  around  solid  bodies  does  not  lead  to  an 
increase  in  the  noise.  This  result  is  also  encouraging  since,  typically,  NS  results  become  noisier  when 
increasing  the  spatial  grid  resolution  near  the  surface. 
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Figure  32  Shear  and  normal  force  for  the  problem  of  gas  flow  around  a  cylinder  for  the  mesh 

refinement  level  of  9  near  the  surface. 


Figure  33  compares  the  calculated  drag  coefficient  with  experimental  data  54.  Also,  results  of  DSMC 
simulations  are  shown  in  Figure  33  for  M=2  and  UFS  simulations  for  M=3.  It  is  seen  that  UFS  results 
agree  well  with  the  experimental  data  and  DSMC  simulations. 
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Figure  33  Drag  coefficient  versus  Kn  number.  Comparison  of  UFS  simulations  with  DSMC 
simulation  and  experiment.  Solid  lines  indicate  the  free  molecular  flow  and  continuum  limits. 

The  problem  of  supersonic  flow  past  cold  bodies  (temperature  less  than  the  stagnation  temperatures) 
is  difficult  since  it  requires  careful  resolution  of  the  Knudsen  layer  around  the  surface.  In  this  layer  the 
density  rises  sharply  and  the  mean  free  path  decreases  strongly.  The  farther  the  temperature  of  the 
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Y  —  1  2 

body  is  from  the  stagnation  temperature  Tstag  =  1  +  ^  M  7^,  the  more  non-equilibrium  the 

distribution  function  is  around  the  solid  body  surface.  For  y  —  5/3 ,  one  obtains  Tstag  =  4  at  M=  3 
and  =  34  at  M=  10. 


We  first  studied  this  problem  for  M=  3  with  a  surface  temperature  Tw  =  1 .  This  case  is  less  demanding 

numerically  than  the  M=  10  case  and  has  been  run  on  a  8-node  cluster.  Figure  34  shows  the  profiles  of 
macroparameters  along  the  centerline.  One  can  see  that  the  temperature  starts  to  drop  (from  a  high 
value  of  35  to  about  5)  at  a  distance  of  about  30  mean  free  paths  (X)  from  the  cylinder  surface  (see 
plot  of  (x-R)A,  in  Figure  34).  This  region  requires  kinetic  description  with  the  Boltzmann  solver. 
Choosing  smaller  kinetic  regions  results  in  non-monotone  profiles  of  macroparameters  and  pressure 
and  heat  transfer  coefficients  (see  Figure  35).  Figure  36  shows  the  computational  grid  and  kinetic 
region  (in  red)  for  the  base  case  for  M=  3. 


Flow  Past  Cylinder:  UFS:  Mach=3,  Kn=0.01,  Tw=1 


Figure  34.  Profiles  of  density ,  velocity  and  temperature  for  3  cases  with  different  numbers  of  kinetic 
cells  for  M  =  3  .  Shown  also  is  the  profile  of  normalized  distance  (x-R)/A  from  the  cylinder  surface. 
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Figure  35.  Profiles  of pressure  and  heat  transfer  coefficients  for  3  cases  with  different  numbers  of 

kinetic  cells  for  M  =  3. 


Figure  36.  Computational  grid  and  kinetic  region  for  the  M=3  base  case. 


The  M=  10  case  is  more  difficult  numerically  since  it  requires  larger  spatial  and  velocity  grids.  This 
case  has  been  run  on  an  HP  cluster  with  larger  number  of  nodes.  For  128  processors,  the  typical  wall 
clock  time  is  14  hrs  and  the  memory  usage  is  710  MB  per  processor.  Figure  37  shows  the  profiles  of 
macroparameters  along  the  centerline  for  the  M=  10  case.  One  can  see  that  the  temperature  starts  to 
drop  (from  a  high  value  of  35  to  about  5)  at  a  distance  of  about  35  mean  free  paths  from  the 
cylinder  surface  (see  plot  of  ( x-R)/X  in  Figure  37).  This  region  requires  kinetic  description  with  the 
Boltzmann  solver.  The  pressure  and  heat  transfer  coefficients  are  shown  in  Figure  38.  The  heat 
transfer  coefficient  is  slightly  noisy  and  can  be  improved  by  further  refining  the  spatial  grid  at  the 
cylinder  surface.  The  current  results  are  obtained  on  a  spatial  grid  in  which  the  cell  size  is  about  2 
mean  free  paths  (X)  near  the  cylinder  body.  The  M=  3  case  shows  that  to  obtain  smooth  heat  fluxes,  the 
cell  size  of  the  order  of  X  is  required.  It  is  worth  noting  that  DSMC  codes  typically  require  the  cell 
size  to  be  smaller  than  local  X. 


Flow  Past  Cylinder:  UFS:  Mach=10,  Kn=0.01,  Tw=2.5 
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Figure  37.  Profiles  of  density ,  velocity  and  temperature  for  M  =  10  case.  Shown  also  is  the  profile  of 
normalized  distance  (x-R)/A  from  cylinder  surface. 


37 


8572/F 


Press,  and  heat  transfer,  Mach=10,  Kn=O.Q1,  Tw=2.5 


Press,  and  heat  transfer,  Mach=lOf  Kn=0.01 ,  Tw=2.5 


Figure  38.  Pressure  and  heat  transfer  coefficients  for  the  M=10  case  (linear  scale  on  the  left  and  log 

scale  on  the  right). 

We  have  compared  the  UFS  results  with  DSMC  results  (SMILE  code)  for  the  case  Mach=10, 
Kn=0.01  case  described  in  55.  The  Euler/BGK  option  was  used  in  UFS  simulations.  The  BGK  solver 
was  used  only  near  the  cylinder  surface  (not  in  the  shock  wave  area).  As  seen  in  Figure  39,  the  UFS 
results  agree  well  with  the  SMILE  results.  As  expected,  the  shock  structure  is  not  resolved  well  with 
the  Euler  solver,  but  the  jumps  of  macro  parameters  are  predicted  well.  The  flow  is  not  fully 
converged  behind  the  cylinder  in  these  simulations. 
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Figure  39  .  Profiles  along  cylinder  surface  (top)  and  along  centerline  (bottom). 
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4.1.2.  Hollow  cylinder  flare  and  biconic  configurations 

The  hollow  cylinder  flare  configuration  has  been  set  up  according  to  Ref.  [56].  The  UFS  calculations 
of  the  simplified  (upper  half)  configuration  are  shown  in  Figure  40.  The  results  are  close 
quantitatively  with  those  presented  in  Ref.  [56]  (see  temperature  and  pressure  distributions  in  Figure 
40).  Further  detailed  study  of  these  cases  is  envisaged. 


Figure  40.  Computational  grid  and  profiles  of  temperature  and  pressure  for  the  hollow  flare  case. 
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Fig,  2  Temperature  contours  TfT.D 


Figure  41.  Pressure  and  temperature  distributions  from  Ref  [56]. 

Figure  42  shows  results  of  preliminary  simulations  of  a  complete  hollow  cylinder  flare  configuration. 
The  computational  grid  inside  the  cylinder  body  is  intentionally  made  coarse.  Detailed  comparison 
with  the  DSMC  results  of  [57]  is  planned. 


Figure  42.  Computational  grid  and  Mach  number  for  the  hollow  flare  case  for  complete  geometry. 


4.2.  Micro  Channels 


The  UFS  has  been  tested  for  simulations  of  internal  flows  in  channels  and  nozzles.  Figure  43  shows 
an  example  of  2D  simulations  of  a  short  channel  for  two  different  Kn  numbers.  The  geometry  and 
flow  conditions  correspond  to  Ref.  58  The  BTE-NS  option  of  the  UFS  solver  with  2nd  order  was  used 
for  these  simulations.  The  boundary  conditions  on  the  left  boundary  are  pin  =  pin  =  1.5 ,  the  boundary 

condition  on  the  right  boundary  is  pout  =  0.5  . 


40 


8572/F 


Figure  43  UFS  simulations  of  2D  micro  channel  Kinetic  and  continuum  domains  (top),  axial  velocity 
(middle),  density  (bottom).  Left:  Kn=01,  density  0.62<r<1.51,  velocity  0.08<v<0.60;  Right: 
Kn=0.01,  density  0.59<  r  <1.5,  velocity  0.1  l<v<0.92 

We  have  performed  a  series  of  calculations  of  gas  flow  in  a  long  channel  for  a  wide  range  of  Knudsen 
numbers.  The  Knudsen  paradox  (the  minimum  of  the  mass  flow  rate)  was  observed  at  Kn~l,  see 
Figure  44. 


Figure  44  The  Knudsen  effect  in  a  relatively  long  channel  (L/d=21)  for  PO/P  1=2  (left)  and PO/P  1=1.5 

(right). 

Different  molecular  models  have  been  compared.  For  example,  the  calculations  by  two  potentials 
(Hard  Spheres  and  Lennard- Jones)  demonstrated  similar  behavior  of  the  main  quantities.  The  flow 
field  structures  are  quite  similar,  and  the  difference  in  a  mass  flow  rate  is  about  5%  (the  flow  rate  is 
smaller  for  the  hard  sphere  molecules). 

4.3.  Nozzle  and  Plume  Flows 

UFS  has  been  sucessfully  used  for  simulations  of  nozzle  and  plume  flows.  Figure  45  shows  results  of 
a  micro  nozzle  simulations  for  two  Kn  numbers.  The  geometry  and  flow  conditions  are  described  in  59 . 
The  BTE-Euler  option  of  the  UFS  with  the  first  order  numerical  scheme  was  used,  the  grid  adaptation 
is  based  on  S  =  log(p)  +  log (u)  . 
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Figure  45  UFS  simulations  of  2D  micro  nozzle  for  Kn=0.01  (left)  and  Kn=0.001  (right)  .Kinetic  and 
continuum  domains  (top),  local  Mach  number  distribution  (middle),  density  (bottom). 

UFS  has  been  demonstrated  for  coupled  simulations  of  nozzle  and  plume  60 .  Figure  46  shows  the  asi- 
symmetric  nozzle  geometry  of  the  NASA  Ames  EAST  facility  and  the  computational  domain  for  the 
plume  simulations.  The  UFS  simulations  were  carried  out  in  the  entire  domain  including  nozzle  and 
plume.  Figure  47  shows  the  contours  of  Mach  number  and  spatial  distributions  of  the  translational  and 
vibrational  temperatures.  The  continuum  solver  was  used  inside  the  nozzle  and  in  the  dense  part  of  the 
plume.  The  low-density  part  of  the  plume  was  calculated  with  the  Boltzmann  solver.  DSMC 
simulations  were  performed  for  the  comparison  in  the  plume  region.  The  results  of  the  UFS  and 
DSMC  simulations  are  compared  in  Fig.  49 
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Figure  46  Nozzle  geometry  and  computational  domain  for  combined  nozzle/plume  simulations 
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Figure  47  Mach  number,  translational  and  vibrational  temperature  contours  for  nozzle  and  plume 
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Figure  48  Calculated  and  experimental  temperatures  along  the  nozzle  axis 


Figure  49  Comparison  ofUFS  and  DSMC  results:  Axial  flow  velocity  along  the  centerline  (left)  and 
radial  distribution  of  the  rotational  temperature  atx=1.15  cm  (on  the  right) 

4.4.  Low  Speed  Flows 

Deterministic  Boltzmann  solver  has  advantages  over  statistical  methods  for  low  speed  flows.  Figure 
50  from  61  illustrates  this  statement:  “Even  after  a  million  time  steps  of  sampling,  the  statistical 
scatter  is  still  so  large  compared  to  the  small  changes  of  flow  quantities  that  meaningful  results  can 
not  be  obtained. . .  In  contrast  to  the  DSMC  method,  the  CFD  method  based  on  Boltzmann  equation  is 
free  from  the  statistical  scatter” 
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Figure  50  Density  contours  for  flow  around  cylinder,  M=0.1,  Kn=0.1.  DSMC  results  on  the  left, 
Boltzmann  results  on  the  right  (after  Morinishi  61  ) 

Figure  51  shows  results  of  UFS  simulations  of  low  speed  flow  around  a  sphere.  Kinetic  and 
continuum  domains  are  shown  in  the  vertical  plane,  the  velocity  contours  are  in  the  horizontal  plane. 


Figure  51  UFS  simulation  flow  around  sphere  at  M=0. 1,  Kn=0. 1 


Figure  52  illustrates  a  low-speed  flow  induced  by  temperature  gradients.  The  nonuniform  boundary 
temperature  distribution  can  induce  flows  in  reactor:  a  significant  flow  (-2-3  m/s)  of  simulated  by 
simple  gas  radicals  flow  over  semiconductor  wafer  expected.  Flow  speed  of  order  of  0.01  thermal 
speed  is  observed.  This  flow  is  absent  according  to  the  traditional  NS  equations  with  slip  boundary 
conditions.  Both  the  direct  Boltzmann  solver  and  the  kinetic  NS  solver  produce  correct  physical 
picture  of  the  flow  shown  in  Figure  52.  The  temperature  T  of  surfaces  goes  from  1.7  to  1  (hot  bottom 
and  cold  top),  bottom  is  symmetry,  top  and  left  boundaries  have  T=  1 .  Note  that  flow  speed  smaller  at 
free  molecular  flow  ( Kn  =  oo)  and  continuum  conditions  (Kn= 0).  Flow  speed  maximum  at  some  Kn  ~ 
0.1. 
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Figure  52.  Temperature  driven  vortex:  temperature  and  velocity  fields  for  three  values  of  Knudsen 
numbers  (Kn=0.01  (\v\max~5E-4),  0.07  (\v\max~0.007),  0.3  (\v\ max-0.005)  from  left  to  right). 
Kinetic  and  continuum  zones  are  shown  in  the  middle  Figure  corner:  blue  -  continuum,  brown  - 

kinetic  zones. 

4.5.  Space  Vehicles  -  Reentry  Problems 

We  have  tested  and  demonstrated  UFS  for  reentry  problems.  Figure  53,  Figure  54  and  Figure  55  show 
results  of  simulations  for  the  Orbital  Reentry  Experiment  (OREX).  Three  types  of  simulations  have 
been  performed:  1)  simulations  using  the  gasdynamic  solver  (2nd  order  Euler  solver)  of  the  flow  at  low 
Knudsen  numbers  (continuum  regime),  see  Figure  53;  2)  simulations  using  the  Boltzmann  solver  at 
high  Knudsen  numbers  (kinetic  regime),  see  Figure  54;  and  3)  simulations  using  UFS  (coupled 
Boltzmann  and  continuum  solvers)  for  moderate  Knudsen  number  (intermediate  regime),  see  Figure 
55. 


Figure  53  Results  of  simulations  of  the  OREX  at  M  =  25  using  the  2nd  order  Euler  solver  in  UFS. 


46 


8572/F 


Figure  54  Results  of  simulations  of  the  OREX  at  M  =  10  and  Kn  =  10  using  the  Boltzmann  solver  in 

UFS. 


Figure  55  Results  of  simulations  of  the  OREX  at  M  =  10  and  Kn  =  0.1  using  a  coupled  kinetic  Euler 
and  Boltzmann  solvers.  Shown  are  the  density  profile  (vertical  plane)  and  computational  grid  and 

kinetic  domain  (in  red)  on  horizontal  plane. 

As  part  of  validation  and  demonstration  of  the  UFS,  we  have  performed  simulations  of  the  Inflatable 
Reentry  Vehicle  Experiment  (IRVE)  using  Euler+Boltzmann  solvers  on  7-processor  Linux  cluster. 
The  flow  conditions  are  for  91  km  altitude:  Kn  =  0.01  and  Mach  =  3.94  62 .  The  flow  is  at  zero  angle  of 
attack.  Figure  56  shows  streamlines,  Mach  number,  and  computational  mesh  (on  the  left)  and  the  gas 
temperature  T  (in  units  of  free  flow  gas  temp)  in  the  vertical  plane,  as  well  as  kinetic  (red)  and 
continuum  (blue)  domains  in  the  horizontal  plane. 
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Figure  56  Streamlines ,  Mach  number,  and  computational  mesh  (on  the  left).  Gas  temperature  in  the 
vertical  plane,  kinetic  (red)  and  continuum  (blue)  domains  in  the  horizontal  plane  (on  the  right). 


4.6.  Unstable  flows 


The  UFS  can  be  a  useful  tool  for  studies  of  unstable  flows,  in  particular,  for 

•  analysis  of  large  scale  structures  by  means  of  kinetic  Euler  or  Navier-Stokes  solvers 

•  kinetic  analysis  of  instabilities  on  the  basis  of  the  Boltzmann  (or  BGK)  solver 

•  multi-scale  description  of  unstable  (turbulent)  flows 

First  results  of  this  direction  have  already  been  obtained.  Figure  57  shows  the  gas  density  distribution 
for  a  subsonic  flow  (M=0.6)  around  a  prism  obtained  with  the  kinetic  Euler  scheme.  A  vortex 
structure  is  formed  in  the  wake  behind  the  prizm,  as  discussed  in 63 


Figure  5  7.  Large  scales  structures  for  gas  flow  around  a  prism  by  the  kinetic  Euler  scheme,  2nd 

order. 

Figure  58  shows  supersonic  flow  around  a  prism  at  an  angle  of  attack  =  3  deg  according  to  the  kinetic 
Navier-Stokes  scheme.  Unstable  structures  appear  in  a  wake  behind  the  prism  at  M=3,  Kn=0. 00001. 
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Figure  58.  Mach  number  (0.009<M<5)  on  the  left ,  gas  temperature  (0.4<T<4.14)  on  the  right. 

A  subsonic  flow  around  a  cylinder  has  been  studied  for  different  regimes  using  the  Boltzmann  (BGK) 
solver.  It  is  known  from  experiments  and  theoretical  investigations  that  for  0<Re<4  the  flow  is 
symmetrical.  In  Figure  59  the  velocity  field  for  a  steady  regime  at  M=0.55,  Kn=0.2  (Re=2.75)  is 
depicted. 


Figure  59.  Velocity  field  for  a  subsonic  flow  at  the  Reynolds  number  Re=2. 75 


For  Reynolds  numbers  in  the  range  0<Re<40  one  can  expect  attached  vortices  behind  a  cylinder.  In 
Figure  60  this  flow  structure  is  shown  for  M=0.55,  Kn=0.015  (Re=37). 


Figure  60.  A  steady  regime  with  a  pair  of  attached  vortices  obtained  with  the  Boltzmann  (BGK) 

solver. 

For  Reynolds  numbers  in  the  range  40<Re<60-100  one  can  expect  the  appearance  of  the  von  Karman 
vortex  street.  Figure  61  shows  the  results  of  the  simulation  for  M=0.6  and  Kn=0.015  (Re=40).  The 
von  Karman  vortex  street  is  formed  during  the  simulations  without  external  flow  perturbations,  in 
distinction  with  the  results  of  64.  Figure  62  shows  the  axial  U-velocity  vs  time  at  different  monitor 
points. 
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Figure  61.  The  longitudinal  velocity  U  and  computational  grid  (on  the  left),  and  stream  lines  (on  the 

right). 


Figure  62.  The  longitudinal  flow  velocity  vs  time  at  different  spatial  points.  Cylinder  of  radius  R=0.1 
is  located  atx=0,  y=0,  the  monitor  points  are  aty=0  and  x=0.2,  0.4,  0.6,  0.8  and  1.0 
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5. 


EXTENSION  TO  MOLECULAR  GASES  AND  REACTIVE  GAS  MIXTURES 


In  this  section  we  describe  UFS  extensions  for  mixtures  of  atomic  ans  molecular  gases 

5.1.  Mixtures  of  atomic  gases 

We  first  implemented  the  BGK-type  collision  integral  described  in  65 .  This  model  is  based  on  the 
following  set  of  equations 
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dt 
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is  a  Maxwellian  distribution.  The  parameters  of  this  distribution  are  defined  by  conservation  of  total 
momentum  and  energy  of  the  mixture  (see  Ref.  65  for  notations): 
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The  frequencies  are  defined  as  va  =  'S^JlpXap  ’  anc^ 
Xafi  =  6%/2 7na1ap{kTalma  +  kTp  t »7/;)12  are  derived  in  [66]. 


the  quantities 


5.1.1.  Temporal  relaxation  in  a  binary  mixture  of  rare  gases  with  disparate  mass 

Initial  conditions  correspond  to  two  Maxwellian  distributions  with  equal  number  densities 
na~np-  0-5 ,  temperatures  Ta=T/3=  1 ,  and  mean  velocities  va  =  —  0.1 .  The  conservation  of 

momentum  implies  that  the  final  equilibrium  velocity  is  given  by: 


v  =  v 


nama  ~npmp 
nama+npmP 


It  is  expected  that  the  approach  to  equilibrium  occur  in  two  stages.  During  the  first  stage 
(Maxwellization),  the  velocity  and  heat  conduction  relax  with  a  characteristic  scale 
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naVp+npVa 
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(41) 
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During  the  second  stage,  the  temperatures  of  the  species  equilibrate  with  a  characteristic  scale 
vT  =  vnip  /  ma  >v  .  This  stage  can  be  described  by  a  two-temperature  hydrodynamics. 

Dimensionless  units  are  used  in  simulations.  The  particle  velocity  is  measured  in  units  of  thermal 
velocity  of  the  first  (light)  species,  vth  ={2kT / Mx^'2 .  The  Knudsen  number  for  a  gas  mixture  is 
defined  as 

Kn  =  Kn0  j nxd2  +  n2d22  yjl  +  mx  /  m  2  /  2  + . . .  j  (42) 

where  dV]  =  (rf.  +  dj)2  / 4  and  =  \l L  .  We  use  common  momentum  space  so  that  components 
with  different  mass  have  different  velocity  range. 

Figure  63,  Figure  64  and  Figure  65  compare  the  mean  particle  velocity,  heat  flux  and  temperature  of 
species  for  HS  and  BGK  models  in  HeNe  and  HeXe  mixtures. 


Figure  63.  Time  dependence  of  the  mean  velocity  for  HeNe  and  HeXe  mixtures. 
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Figure  64.  Time  dependence  of  species  temperature  for  HeNe  and  HeXe  mixtures. 


Figure  65.  Time  dependence  of  heat  fluxes  for  HeNe  and  HeXe  mixtures. 

It  is  seen  from  Figure  63,  Figure  64  and  Figure  65  that  relaxation  processes  in  gas  mixtures  with 
disparate  mass  occur  in  two  stages  with  vastly  different  time  scales.  The  first  (short)  stage 
corresponds  to  Maxwellization  of  the  distribution  function.  The  heat  flux  vanishes  at  the  end  of  this 
stage.  The  second  (long)  stage  corresponds  to  relaxation  of  temperatures  to  a  common  equilibrium 
value.  This  stage  can  be  described  by  a  two-temperature  hydrodynamic  model. 
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5.1.2.  Spatial  Relaxation  for  a  Gas  Mixture 


The  nonuniform  spatial  relaxation  problem  67  is  an  analog  of  the  temporal  relaxation  problem  and 
describes  the  relaxation  processes  in  space.  At  the  boundary  of  a  half-infinite  region  (x=0)  the 
nonequilibrium  distribution  function  is  specified.  The  nonequilibrium  distribution  is  relaxed  to 
equilibrium  downstream  at  x>0.  The  Mach  number  is  assumed  to  be  sufficiently  large  thus  the 
negative  flux  of  the  back- scattered  molecules  can  be  neglected.  The  interesting  properties  of  the 
spatial  behaviour  of  macroparameters  have  been  observed  such  as  anomalous  heat  transfer.  For 
mixtures  of  gases  the  nonequilibrium  structure  can  be  rather  complicated.  This  problem  can  serve  a 
basic  model  for  the  description  of  kinetic  processes  in  open  systems. 

We  first  considered  a  binary  mixture  of  simple  gases.  The  anomalous  heat  transfer  noted  above  for  a 
monatomic  gas  is  also  observed  in  a  binary  gas  mixture.  Namely,  the  gradient  of  the  total  temperature 
can  have  the  same  sign  as  the  total  heat  flux.  This  peculiarity  of  the  nonequilibrium  heat  convection 
denotes  that  the  nonequilibrium  boundary  distribution  can  heat  up  (or  cool  down)  the  region 
downstream,  so  the  temperature  can  increase  (decrease)  if  the  heat  flux  is  positive  (negative) 
respectively.  Figure  66  shows  the  spatial  distribution  of  the  total  temperature  and  the  heat  flux.  The 
spatial  distance  is  given  in  units  of  the  mean  free  path.  The  BGK-type  model  is  used.  Solutions  of  the 
Boltzmann  equation  for  hard  sphere  molecular  model  (HS)  demonstrate  a  similar  behaviour  (the 
relaxation  region  is  larger  than  that  for  the  BGK  model).  In  our  simulations,  the  boundary  densities, 
masses  and  diameters  of  the  gas  components  are  equal  to  the  ratios  for  oxygen  and  nitrogen  in  air.  The 
densities  are  1  and  3.71,  the  masses  are  1.14  and  1  and  the  diameters  1  and  1.07,  respectively  (here 
molecules  of  oxygen  and  nitrogen  is  considered  as  simple  and  unstructured).  The  boundary  conditions 
for  each  component  correspond  to  different  Maxwellian  distributions.  It  is  interesting  that  for  the  1-st 
component  the  temperature  decreases  monotonically  downstream  and  the  heat  flux  in  negative  (an 
absolute  value  of  the  heat  flux  tends  to  zero  downstream),  for  the  2nd  component  the  temperature 
increases  downstream  and  the  heat  flux  is  positive.  It  is  important  that  the  total  temperature  gradient 
and  the  heat  flux  have  the  same  signs,  namely,  this  gradient  and  the  flux  are  both  negative. 


Figure  66.  Spatial  distribution  of  the  total  temperature  (on  the  left)  and  the  total  heat  flux  (on  the 
right)  for  the  nonuniform  relaxation  problem  in  a  binary  gas  mixture 

The  structure  of  the  ralaxation  zone  for  a  multiple  gas  mixture  can  be  more  complex;  in  particular, 
one  can  change  the  spatial  distributions  of  density  by  changing  the  boundary  conditions  for  the  gas 
components.  Figure  67  shows  the  spatial  distributions  of  the  species  density  and  the  total  density  for  a 
mixture  of  3  gases.  For  a  one-component  gas,  the  density  can  only  decrease,  but  for  a  mixture,  the 
total  density  increases  with  distance  (see  Figure  67). 
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This  spatial  relaxation  problem  in  gas  mixtures  has  also  been  studied  for  an  unsteady  case  and  for  a 
2D  steady  setting. 


Figure  67.  Spatial  distribution  of  specie  densities  for  a  nonuniform  relaxation  problem  in  a  3- 

component  mixture 


5.1.3.  Gas  acceleration  by  optical  forces 

Laser  induced  optical  lattices  have  recently  attracted  considerable  interest  for  several  applications 
including  contact-free  diagnostics  68.  The  direct  Boltzmann  solver  has  advantages  over  statistical 
methods  for  the  analysis  of  low  amplitude  density  perturbations  induced  by  laser  fields  We  have 
solved  a  set  of  kinetic  equations  in  the  form 


cv 


(43) 


where  F(x)  =  F0  sin(2^v/Z)  .  Computational  domain  is  from  0  to  L,  and  cyclic  boundary 

conditions  are  set  up  at  x=0  and  x=L.  We  have  studied  the  BGK  and  HS  models  of  inter-molecular 
interactions. 


Figure  68  and  Figure  69  show  the  density  modulation  of  different  species  in  HeNe  mixture  induced  by 
optical  lattices.  The  observed  resonance  corresponds  to  sound  wave  propagation  in  the  binary  mixture 
with  the  sound  speed  determined  by  gas  composition.  The  heavy  and  light  components  respond 
differently  to  laser  forces.  Figure  68  (right  part)  illustrates  the  effect  of  gas  composition  on  density 
perturbations,  reflecting  the  change  in  sound  speed  in  gas  mixtures  of  different  composition..  Figure 
69  shows  the  effect  of  Kn  number  on  the  density  perturbations 
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HeNe:  0.5:0.5,  Kn=0.01 


HeNe:  BGK,  Kn=0.01 


Figure  68.  Density  modulation  of  different  species:  comparison  of  BGK  and  HS  models  (left).  Total 
density  according  to  BGK  model  for  different  fraction  of  species  (right). 

HeNe,  0.5:0.5,  BGK 


Figure  69.  Density  modulation  of  different  species  for  two  Kn  numbers  via  BGK  model 

These  results  illustrate  the  potential  of  UFS  for  simulations  of  kinetic  effects  induced  by  optical  forces 
in  gas  mixtures. 

5.1.4.  Shock  wave  structure  in  binary  gas  mixture 

Figure  70  shows  results  of  simulations  of  the  shock  wave  structure  for  M  =  3,  the  mass  ratio  m2/ml  = 
3,  and  density  ratio  n2/nl=0.\.  The  boundary  conditions  are  taken  from  Ref.  69 .  The  velocity  mesh 
used  in  these  simulations  had  32*32*16  nodes. 
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Shock  Tube:  MultiSpecies  BGK,  m2/m1=3,  Species  1 
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Shock  Tube:  MultiSpecies  BGK,  m2/m1=3,  Species  2 
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Figure  70  Shock  structure  for  a  two-component  gas  mixture  with  a  mass  ratio  of  3. 

Two  types  of  velocity  meshes  have  been  tested,  the  one  with  same  velocity  mesh  for  different 
components,  and  the  other  with  different  velocity  mesh  (cell  numbers  and  velocity  range)  for 
components  with  different  mass.  The  second  type  of  mesh  turned  out  to  be  preferable  since  it  allowed 
us  to  capture  the  key  behavior  of  the  distribution  function  with  smaller  number  of  mesh  points. 
Criteria  for  selection  of  optimal  mesh  in  velocity  space  were  analyzed. 

Furthermore,  the  collision  integral  based  on  HS  model  for  multi-component  mixtures  70  was 
implemented  using  momentum  space  rather  than  velocity  space.  We  have  carried  out  calculations  for 
the  shock  wave  structure  in  a  mixture  of  two  gases  with  mass  ratio  %,  the  upstream  density  of  heavy 
component  0.9  and  M=2.  The  HS  model  is  used  with  equal  molecular  diameters.  Figure  71  shows  the 
velocity  distribution  functions  (averaged  over  y-  and  z-directions)  of  both  species  on  the  common 
momentum  space  and  at  different  locations  in  the  shock  wave.  Figure  72  shows  the  distribution  of 
normalized  densities  [(nfx)  —  nt_)  /(ni+  -  nt_  )  ]  and  parallel  and  perpendicular  temperatures  of  both 

gas  species.  The  results  are  in  close  agreement  with  70  One  can  see  that  the  parameters  of  the  heavy 
component  react  to  shock  wave  with  a  delay  compared  to  those  of  the  light  component.  The 
temperature  of  the  heavy  component  “overshoots,”  which  is  also  a  well  known  phenomenon  70.  We 
also  note  that  with  increasing  the  mass  ratio,  the  computational  time  increases  greatly.  Using  different 
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velocity  mesh  for  different  components  and  simplifying  collision  integral  for  the  mixtures  with 
disparate  mass  71  should  increase  the  efficiency  of  simulations. 


Figure  71  Velocity  distributions  of  light  (top)  and  heavy  (bottom)  species  at  different  points  of  the 
shock  wave  for  M=2  and  mass  ratio  The  species  momentum  space  px  is  used  as  x-axis. 


SW/HS/M=2/m1:m2=4:1 


SW/HS/M=2/m1  :nn  2=4:1 


Figure  72  Profiles  of  normalized  density,  velocity,  temperature  and  parallel  and  perpendicular 
temperatures  for  the  shock  wave  in  a  binary  gas  mixture  for  M=2  and  m1/m2  = 
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Figure  73  shows  the  results  of  simulations  using  the  multi-species  BGK  Boltzmann  solver  for  multi- 
component  gas  flow  around  a  cylinder  at  Mach  =  3  and  Kn  =  0.001  with  3  species  of  the  following 
masses:  mi  =  1,  m2  =  1.5  and  m3  =  2.  The  collision  diameters  du  d2  and  d3  assumed  to  be  the  same  and 
equal  to  1 . 


Figure  73  Multi-species  BGK  Boltzmann  solver:  spatial  distributions  of  total  gas  density  (top, 
0.24<rho<8.97)  and  temperature  (bottom,  1<T<5.3)  for  M  =  3,  Kn  =  0.001,  3  species  of  masses:  mj 
=  1,  m2  =  1.5  and  m3  =  2,  with  equal  collision  diameters. 

5.1.5.  Multi-species  Euler  solver 

The  multi-species  continuum  solver  has  been  implemented  in  UFS.  The  code  enables  user  to 
automatically  create  any  number  of  variables  for  solving  continuum  equations  for  each  species  and 
specify  initial  conditions,  boundary  conditions  at  the  boundaries  of  computational  domain  and  at  solid 
body  surfaces.  The  implementation  is  based  on  the  multi-species  BGK  operator  with  an  assumption  of 
instantaneous  equilibration  of  all  mixture  components. 

Figure  74  shows  results  of  simulations  using  the  multi-species  kinetic  Euler  solver  for  a  flow  around  a 
cylinder  at  M  =  3,  with  3  species  of  the  masses:  mi  =  1,  m2  =  1.5  and  m3  =  2,  with  equal  collision 
diameters  di?  d2  and  d3  =1. 
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Figure  74  Multi-species  kinetic  Euler  solver:  Spatial  distributions  of  total  gas  density  (top, 
0.22<rho<10.5)  and  temperature  (bottom,  1<T<4.79)  for  M  =  3,  Kn  =  0.001,  3  species  of  masses: 
mi  =  1,  m2  =  1.5  and  m3  =  2,  with  equal  collision  diameters. 

The  coupling  algorithm  between  the  multi-species  Euler  and  Boltzmann  solvers  has  been  developed 
for  arbitrary  gas  mixtures. 

5.2.  Rotationally  Excited  Molecules 

5.2.1.  Wang-Chang-Uhlenbeck  (WCU)  solver 

We  have  developed  a  Wang-Chang-Uhlenbeck  (WCU)  solver  for  molecular  gases  with  rotational 
degrees  of  freedom  following  72 .  The  WCU  solver  was  tested  for  the  shock  wave  (SW)  structure  in 
Nitrogen  for  a  wide  range  of  Mach  numbers  and  compared  with  experimental  cases  of  Alsmeyer  (for 
M=  1.53,  1.7,  2,  2.4,  3.2,  3.8,  6.1,  8.4,  10)  and  Robben  and  Talbot  experiment 73  for  M=1 . 

Figure  75  compares  experimental  data  by  Alsmeyer  with  simulation  results  for  M=1.53  and  M=l.  7. 
The  computations  were  performed  using  30  rotational  levels  for  the  first  case  and  32  levels  for  the 
second  case.  The  number  of  levels  was  selected  based  on  the  temperature  range  of  the  considered 
problem.  For  small  deviations  from  the  room  temperature,  using  24  levels  for  Nitrogen  was  found  to 
be  sufficient.  The  typical  CPU  time  was  about  35  hours  in  both  cases  on  a  AMD64  3000  processor. 


Figure  75  Shock  wave  structure  in  Nitrogen  for  M=  1.5 3  (left)  and  M=l.  7  (right) 

Figure  76  shows  rotational  spectrum  for  44  levels  at  several  points  of  the  wave  front: 
v  =  -oo?  x  =  X  ,  -  2/1,  x  =  vc  -  A,  x  =  xc ,  v  =  xc  +  A,  x  =  oo  .  Here  xc  denotes  the  SW  center  defined 

as  the  point  where  the  reduced  gas  density  is  equal  to  Vi.  On  the  x-axes  is  the  number  of  the  rotational 
level,  on  the  y-axes  is  the  population  of  the  nth  rotational  level.  It  is  seen  that  some  peculiarity  of  the 
spectrum  around  7-8  levels  is  observed  at  the  second,  third  and  fourth  spatial  points. 
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Figure  76  Rotational  spectrum  at  different  points  of  the  shock  wave  atM=3.2. 

Figure  77  shows  distributions  of  gas  density,  translational  and  rotational  temperatures  obtained  in  our 
simulations  for  M=  12.9.  Our  results  are  in  good  agreement  with  the  experimental  data  by  Robben  and 
Talbot  and  computations  of  K.  Koura  by  DSMC  method  74. 


Figure  78  shows  rotational  spectrum  for  25  levels  at  several  points  along  the  wave  front.  The  center  of 
SW  is  located  at  X=0.  On  the  x-axes  is  the  rotational  level  number,  on  the  y-axes  is  the  population  of 
the  rotational  levels.  It  is  clearly  seen  that  the  rotational  equilibrium  inside  the  SW  doesn’t  exists  for 
such  a  high  Mach  number. 
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Figure  78  Rotational  spectra  at  different  points  along  the  shock  wave  in  Nitrogen  at  M=12.9. 

The  WCU  model  for  rotationally  excited  molecules  is  rather  expensive  computationally  because  a 
separate  kinetic  equation  is  solved  for  each  rotationally  excited  level  of  the  molecule.  For  practical 
applications,  we  have  implemented  a  simplified  Rykov  model,  which  was  demonstrated  in  a  number 
of  papers  (see  75  for  further  references) 

5.2.2.  Rykov  Model 

Rykov’ s  model  (R-model)  for  molecular  gases  with  rotational  degrees  of  freedom  introduces  two 
functions  f0  =  J  fde  and  f  =  J  efde ,  where  /  is  the  distribution  function  depending  on  rotational 
energy  e .  Then  one  obtains  two  equations 


(44) 


(45) 


The  right-hand  side  of  each  equation  corresponds  to  a  model  collision  integral.  The  macroscopic 
parameters  (corresponding  to  rotational  and  translational  movement  of  molecules)  are  obtained  by 
integration  of  /0  and  f  over  velocity  space  (details  can  be  found  in  [75]).  The  R-model  collision 

integral  has  been  implemented  for  3D  and  2D  cases  (2D  case  takes  into  account  the  symmetry  of  the 
velocity  space).  The  conservation  of  mass,  momentum  and  total  energy  is  enforced  by  introducing  a 
procedure  analogous  to  conservative  correction  used  in  the  BGK  solver. 

The  R-model  was  tested  for  rotational  relaxation  in  nitrogen.  Figure  79  shows  the  time  dependence  of 
the  macroparameters  during  the  relaxation  of  rotationally  excited  nitrogen  molecules. 
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Nitrogen  relaxation  with  R-model 


Figure  79.  Time  dependence  of  the  macroparameters  during  the  relaxation  of  rotationally  excited 

nitrogen 


5.3.  Vibrationally  Excited  Molecules 

We  have  extended  the  WCU  solver  for  vibrationally  excited  molecules.  The  results  for  vibrational 
relaxation  are  shown  in  Figure  80  for  the  initial  conditions:  rv=0.1  and  Tt=\.  It  is  seen  that  the 
translational  and  vibration  temperatures  relax  to  the  same  values  with  the  characteristic  relaxation 
time  depending  on  the  V-T  cross  sections.  The  total  V-T  cross  section  is  higher  then  the  partial  cross 
sections  because  of  a  number  of  permitted  transitions.  The  problem  is  solved  for  12  vibration  levels, 
rotational  levels  are  ignored. 


Figure  80  V-T  relaxation  for  P?  =  0.001,  (kj  ^  /,  j)  and  for  P?  =  0.0001,  (kj  ^  /,  j)  (right) 

5.3.1.  Three-temperature  BGK  model 

In  order  to  account  for  both  the  rotational  and  vibrational  energies  in  the  Boltzmann  solver,  we  have 
developed  a  3-temperature  (or  3T)  BGK  model.  This  model  is  based  on  the  same  ideology  as  the 
previously  implemented  Rykov  model  for  the  rotational  energies. 

In  the  presence  of  internal  degrees  of  freedom,  the  kinetic  equation  can  be  written  in  the  form: 
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(46) 


+  v  .(so  =  r'+r 

dt 

where  Iel  represents  the  elastic  collision  integral  and  Iin  denotes  the  inelastic  collision  integral.  For 
the  treatment  of  inelastic  collisions,  we  use  a  simplified  relaxation  model.  The  elastic  collision 
integral  is  written  in  the  BGK  form: 


where 


{nTf\nTrfr,\nTvfJ1 


I  -V  if m-f)> 

fm=f:(n,u,Tt)f:(Tr)f:(Tv)= 

exp[-(£  -  uf  /  Tt  -  Kr&  K2T)  -  K&  /(2T)\ . 


(47) 


(48) 


Here  Tt  ,  Tr  ,  Tv  denote  the  translational,  rotational,  and  vibrational  temperatures  defined  from  the 
relations: 


nEt  =\.5nTt  +  nu2  =  J  fmf2dfdfdf 

(49) 

nEr  =  0.5nKrTr  =  J  fj2dfdfdf 

(50) 

nEv  =  0.5nKvTv  =  J  fj2dfdfdf , 

(51) 

Here  Kr  and  Kv{T)~ 


29 


—  are  the  numbers  of  rotational  and  vibrational  degrees  of 
rexp(0/r-i) 

freedom  ( 9  denotes  the  characteristic  vibrational  temperature). 

The  inelastic  collision  integral  is  presented  in  the  form: 

rn=v‘r\f:-f)+vtr(f:-f), 

where 

- exp[-(£  -uf/Teq-  Krg  1(2 Teq )  -  Kftf  /(2Teq )] , 


(52) 


(nTeq)  r 


3/2+K/2+KJ2 


(53) 


€  =f>,u,rqr)f:(Teqr)fm(Tv)  = 


-exp[-(^  -  uf  /  Vqr  -  Kr&  /(2  Teqr)  -  Kvg  1(2  Tv)] , 


(7rTeqrfn+K'n(7rTv)KJ2 

Here  Teq  is  the  equilibrium  temperature  of  translational,  vibrational  and  rotational  degrees  of 
3  t  _| ~  j1  -\-  K  T 

freedom,  Teq  =  — L - — - 5LJL ,  and  Teqr  is  the  equilibrium  temperature  of  rotational  and 

3 +k+k:9 


(54) 
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translational  degrees  of  freedom,  Teqr  =  — - - l-L- . 

3  +  Kr 

For  the  numerical  solution,  Eq.  (46)  is  reduced  to  a  set  of  equations  for  three  functions 

j\=\f£dZrdZv, 
f2=\  fgd&dZ, . 


(55) 

(56) 

(57) 


Such  a  reduction  is  possible  because  the  convection  part  of  the  kinetic  equation  does  not  contain 
derivatives  with  respect  to  vibrational  and  rotational  energies  of  the  molecules.  After  integration  of 

Eq.  (46)  with  weights  (1,  <^2,  <^2)  ,  we  obtain: 


dt 


=  vel(-fo  +  fl(Tt))  +  v,rv(-fo  +  X(T«))  +  v\-U  +  C(Teqr)) 


(58) 


=  ve/(-./;  +  EJJTJ)  +  v'n'(-/1  +  Eerqf^(Teq))  +  v\-fx  +  Eerqrf:(Teqr)) 
dt 


(59) 


dt 


ve'(-f2  +  EX<Tt ))  +  v"t(-/2  +  Eeqfl(Teq))  +  vtr(-f2  +  EX(rqr)) 


(60) 


The  total  (elastic  and  inelastic)  collision  frequency  v  is  defined  as  v  =  vel  +  vtrv  +  vtr .  By 
introducing  vtrv  —v / Zv ,  =  v / Zr ,  we  obtain  ve/  =  v(l  —  1/ZV  —  1/Zr),  where  Zr  and  Zv 

define  the  number  of  rotational  and  vibrational  collisions  per  one  elastic  collision  (which  can  be  any 
functions  of  macroparameters  such  as  temperature). 

The  set  of  3  equations  (58-60)  resembles  the  Landau-Teller  model  for  the  relaxation  of  rotational  and 
vibrational  temperatures.  In  the  UFS  code,  we  have  implemented  two  options  for  defining  Z-numbers. 
One  can  specify  them  as  constants  or  use  experimental  rates  (such  as  the  Millikan  rates). 

We  have  further  implemented  and  tested  this  model  in  the  kinetic  3T-Continuum  solver.  The  coupling 
of  the  3T-BGK  and  3T-Continuum  solvers  has  been  tested  and  demonstrated  for  the  ID  Shock  Wave 
(SW)  and  2D  plume  problems. 

5.3.2.  Shock  Wave  in  Nitrogen 

The  developed  3T-UFS  has  been  benchmarked  for  ID  SW  simulations  using  the  results  presented  in 
76.  In  that  work,  Mach  =  5  SW  has  been  simulated  for  nitrogen  gas  taking  into  account  the  rotational 
and  vibrational  energies  using  the  BGK-NS  solver  and  the  DSMC  solver.  The  values  of  Zr  =  3  and  Zv 
=  100  were  used  in  that  paper  and  in  our  simulations.  We  first  carry  out  simulations  using  the  3T- 
Continuum  Solver,  the  results  for  which  are  presented  in  Figure  81.  One  can  see  that  the  3T- 
Continuum  Solver  result  give  sharp  jumps  of  the  macroparameters  as  the  SW  location  and  slow 
evolution  of  Trot  and  Tvib  towards  equilibrium  with  Trot  reaching  equilibrium  quicker  than  Tvib.  The 
results  of  the  3T-Continuum  Solver  are  close  to  those  of  the  BGK-NS  presented  in  [76]  considering 
that  the  Sutherland  law  for  viscosity  has  been  using  in  that  paper  with  unclear  parameters. 
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SW:  3T-Euler  Model,  N2,  Mach  =  5,  Zr  =  3,  Zv  =  100 
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Figure  81.  ID  shock-wave  profiles  of  macroparameters  obtained  using  the  developed  3T-Continuum 

solver. 

We  next  carry  out  simulations  for  the  same  conditions  using  the  developed  3T-BGK  model  in  the 
Boltzmann  solver.  The  results  are  presented  in  Fig.  82.  One  can  see  that  now  there  is  significant 
spread  of  SW  parameters  into  the  upstream  direction  in  particular  of  Tvib.  This  spread  of 
macroparameters  is  close  to  that  obtained  using  the  DSMC  simulations  in  [76].  Also,  the  parallel  and 
perpendicular  temperatures  (their  difference  demonstrates  departure  from  translational  equilibrium) 
are  close  to  those  presented  in  [76]. 
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SW:  3T-BGK  Model,  N2,  Mach  =  5,  Zr  =  3,  Zv  =  1 00 


Figure  82.  ID  shock-wave  profiles  of  macroparameters  and  parallel  and  perpendicular  temperatures 
obtained  using  the  Boltzmann  Solver  with  the  3T-BGK  model. 


We  finally  demonstrate  the  UFS  capabilities  for  this  problem  by  using  coupled  solution  of  the  3T- 
BGK  solver  with  the  3T-Continuum  solver.  The  results  are  presented  in  Figure  83.  The  Boltzmann 
solver  was  used  only  in  the  upstream  region  and  in  the  region  around  the  SW  (see  kinetic  flag  in 
Figure  83).  The  3T-Continuum  Solver  was  used  in  the  extended  downstream  region  where  slow 
evolution  of  Tvib  takes  place.  One  can  see  that  the  results  are  very  close  to  those  obtained  using  the 
3T-BGK  Boltzmann  Solver  with  the  UFS  run  requiring  almost  4  times  less  memory  and  having  much 
faster  convergence. 


SW:  3T-UFS  Model,  N2,  Mach  =  5,  Zr  =  3,  Zv  =  100 
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SW:  3T-UFS  Model,  N2,  Mach  =  5,  Zr  =  3,  Zv  =  100 


Figure  83.  ID  shock-wave  profiles  of  macroparameters,  kinetic  flag,  and  parallel  and  perpendicular 
temperatures  obtained  using  UFS  with  coupling  of  the  3T-BGK  and  3T-Continuum  solvers. 


5.3.3.  3  Temperature  Model  for  a  mixture  of  molecular  and  atomic  gases 


When  more  than  one  species  present  in  a  gas,  one  needs  to  take  into  account  momentum  and  energy 
exchanges  between  them.  If  there  are  chemical  reactions,  the  model  needs  to  take  into  account  also 
changes  in  mass  density  of  each  species.  The  momentum  and  energy  exchanges  between  molecular 
and  atomic  species  in  a  gas  are  taken  into  account  according  to  the  Morse  and  Oguchi  model 
described  in  Ref.  [77].  The  collisional  integral  for  a  mixture  of  atomic  A  and  molecular  A2  species  is 
written  in  the  following  form 


£>*/,,=  Z  <»K,-4)+<» 

b=a,a2 

DU  A  =  Vaa  -fA)  +  <  (M^2  - 


(KB-fA2)+<B(K2S-fA2) 


(61) 

(62) 


where  the  Ms  denote  the  Maxwellian  distribution  functions  calculated  at  corresponding  equilibrium 
temperatures  and  each  term  represents  different  type  of  collisional  process. 


The  validation  of  the  developed  3T-BGK  model  has  been  carried  out  by  simulating  flow  past  cylinder 
for  N2  gas  for  different  Mach  and  Knudsen  numbers.  The  results  presented  in  Figure  84  and  Figure  85 
are  for  Mach  =  3  and  Kn  =  0.05.  One  can  see  in  Figure  85  that  the  vibrational  temperature  starts  to 
grow  well  behind  the  shockwave  and  that  the  3  temperatures  equalize  downstream  far  from  the 
cylinder. 
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Figure  84.  Flow  past  cylinder,  Mach  =  3,  Kn  =  0. 05.  Shown  are  the  computational  mesh  and  the 
kinetic  flag  (in  red  color),  the  Mach  number  profiles  and  the  profiles  of  translational,  rotational  and 

vibrational  temperatures. 
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Flow  Past  Cylinder:  3T-UFS:  N2,  Mach  =  3,  Kn  =  0.05 


x,  distance  along  centerline 

Figure  85.  Flow  past  cylinder  for  Mach  =  3  and  Kn  =  0.05.  Shown  are  the  profdes  along  the 
centerline  of  the  kinetic  flag,  density,  velocity,  translational,  rotational  and  vibrational  temperatures. 


The  3T-BGK  model  has  also  been  tested  for  the  Shock  Wave  (SW)  problem.  Figure  86  shows  the 
results  for  the  SW  problem  at  Mach=2  for  a  mixture  of  a  molecular  and  an  atomic  gases. 
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Figure  86.  Profiles  of  density,  velocity,  translational  temperature,  rotational  and  vibrational 
temperature  for  the  molecule  for  a  ID  SW  in  a  mixture  of  atomic-molecular  gases. 
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5.3.4.  Hypersonic  Nitrogen  flow  past  a  blunt  body 

The  low  enthalpy  case  for  a  flow  past  an  axisymetric  blunt  body  at  Mach=16  has  been  investigated. 
The  case  is  described  in  Ref.  78.  For  these  conditions,  the  UFS-3T  model  was  used  which  takes  into 
account  the  internal  energies  of  molecules.  Model  N2  molecule  was  used  for  simulations,  which 
approximates  the  airflow  experimental  conditions.  The  geometry  of  the  body  and  2D  distrubution  of 
flow  parameters  is  shown  in  Fig  87.  The  flow  parameters  along  the  streamline  are  shown  in  Fig.  88. 
The  shock  wave  position  agrees  well  with  NS  simulations  presented  in  Ref.  [78].  Also,  the  pressure 
distribution  obtained  using  UFS  agrees  well  with  the  NS  results  and  experimental  results  from  Ref. 
[78]  (see  also  Figure  89). 


body. 
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Figure  88.  Profiles  of  macroparameters  along  stagnation  streamline  for  Mach= 16  flow  past  a  blunt 

body. 
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Blunt  Body:  Mach  =  16.34 


Figure  89.  Comparison  of  static  pressure  results  obtained  with  UFS  and  experimental  data  presented 

in  Ref  [78]. 

We  have  also  carried  out  simulations  for  a  plain,  2D  geometry.  The  results  in  Figure  90  show  that  the 
shock  wave  stands  much  farther  from  the  cylinder  surface  compared  to  the  axi-symmetric  case. 


Flow  Past  Cylinder:  UFS-3T:  Mach  =  16.34,  Kn  =  1.67E-04 


r)/Rn,  distance  along  stagnation  streamline 

Figure  90.  Profiles  of  macroparameters  along  stagnation  streamline  for  Mach= 16  flow  past  circular 

cylinder  for  a  2D,  plain  case. 

Details  regarding  UFS  evaluation  for  computing  heat  flux  in  hypersonic  flows  can  be  found  in  79 .  It 
was  concluded  that  using  kinetic  solvers  (Boltzmann  solver  and  CFD  gas  kinetic  schemes)  in 
combination  with  Cartesian  mesh  offers  notable  improvements  of  the  heat  flux  calculations  compared 
to  the  traditional  NS  solvers  and  comparable  to  DSMC  accuracy.  Boltzmann  solver  gives  smooth  heat 
fluxes  for  spatial  grids  with  cell  sizes  of  the  order  of  local  mean  free  path  or  larger,  provided  that 
spatial  grids  with  gradual  stretching  are  used  near  the  surfaces  of  the  body. 

The  optimization  with  respect  to  the  size  of  the  kinetic  domain  (on  a  given  spatial  grid)  has  also  been 
studied  in  [79].  It  was  observed  that  for  a  small  size  of  the  kinetic  domain,  unphysical  (e.g.,  non¬ 
monotone)  profiles  can  be  obtained  near  the  body.  For  an  optimal  size  of  the  kinetic  domain,  reliable 
heat  fluxes  are  obtained  with  minimal  computational  cost.  It  was  also  observed  that  using  the  kinetic 
Navier-Stokes  solver  instead  of  the  kinetic  Euler  solver  allows  reducing  the  size  of  kinetic  domain. 
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5.4.  Chemical  Reactions 


A  chemistry  module  has  been  added  to  the  kinetic  Euler  solver  for  reacting  gas  mixtures.  In  this 
particular  module,  17  reaction  steps  have  been  implemented  between  major  air  species  O,  02,  N,  N2 
and  NO.  Figure  91  shows  results  of  a  relaxation  problem  with  the  following  set  of  chemical  reactions: 

02+M< - >20  +  M 

N2+M  < - >2  N  +  M  (63) 

NO  +  M< - >N  +  0  +  M 

where  M(02,0,N2,N ,NO) .  The  initial  densities  /?0of  species  are  (  in  kg/m3): 
02=l.E-3;0  =  2.E-5;N2=2.E-3,N  =  2.E-4,N0  =  2.E-3 


0.5  r 


o 

Cl 


NO 

N2 


Figure  91  Time  evolution  of  gas  species  due  to  chemical  reactions. 


We  have  extended  the  multi-component  kinetic  Euler  solver  for  reactive  gas  mixtures  following  the 
work  80  taking  into  account  changes  in  the  internal  energy.  The  model  assumes  that  the  translational 
and  rotational  temperatures  are  the  same  and  that  each  molecular  species  has  its  own  vibrational 
temperature  Tv.  Figure  92  shows  results  of  2D  simulations  for  supersonic  air  flow  around  a  cylinder  at 
M=  2.  The  incoming  gas  density  is  p  =  0.1  kg/m3.  It  is  known  that  due  to  slow  V-T  relaxation, 
vibrational  temperature  lags  behind  the  translational  temperature  and  in  our  case  remains  high  behind 
the  cylinder.  Figure  92  shows  distributions  of  species  density,  w-velocity,  translational  and  vibrational 
temperature  of  different  molecules  along  the  central  streamline.  It  is  seen  that  the  vibrational 
temperature  of  different  molecules  is  different  at  low  gas  density. 
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M=2/AirChemistry/pref  =  0.1  kg/m3 


M=2/Air  Chemistry/pref  =  0.1  kg/m3 
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Figure  92  Distributions  of  species  density,  u-velocity,  translational  and  vibrational  temperature  of 
different  molecules  along  the  central  streamline. 


In  order  to  validate  UFS  for  ID  SW  structure,  we  have  carried  out  simulations  similar  to  those 
presented  by  Surzhikov  81  for  ID  SW  structure  in  reacting  air.  Simulations  are  performed  only  for  the 
region  behind  the  SW  since  it  is  difficult  to  impose  correct  boundary  conditions  at  the  SW  boundaries 
for  a  reacting  mixture  taking  into  account  the  vibrational  energies.  The  results  in  Ref.  [81]  are  for 
high-speed  flows  of  11.36  km/s  when  the  temperature  in  the  SW  jumps  up  to  60000  K.  We  have 
carried  out  simulations  of  SW  structure  behind  the  SW  front  for  M  =  15  and  temperature  jump  of 
about  10,000  K.  The  results  are  presented  in  Figure  93.  One  can  see  that  all  3  vibrational  temperatures 
Tv  equilibrate  within  a  distance  of  about  10  cm.  The  species  densities  also  change  over  the  same 
distance.  Simulations  for  higher  temperatures  (such  as  those  in  Ref.  [81])  will  require  a  more  robust 
chemistry  module,  which  will  account  for  very  high  values  of  the  reaction  rates  at  large  temperatures. 
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SW  with  Chem/M=1 5/0.003  kg/m3 


SW  with  Chem/M=1 5/0.003  kg/m3 
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Figure  93  Distributions  of  species  density,  u-velocity,  translational  and  vibrational  temperature  of 
different  molecules  for  ID  SW  structure  behind  the  shock  at  Mach  =15  . 


Development  of  general-purpose  chemistry  in  the  UFS  framework  is  planned  in  the  future.  We  have 
evaluated  the  CANTERA  code  82  which  has  established  itself  as  a  reliable  and  user  friendly  chemistry 
modeling  tool  with  extensive  database  for  various  gas  mixtures.  We  have  developed  an  interface 
between  a  UFS-like  C  code  and  the  CANTERA  code  and  demonstrated  feasibility  of  coupling  UFS 
and  CANTERA  codes  in  the  future. 
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6. 


CONCLUSION 


During  this  Phase  II  SBIR  Project,  we  have  developed  a  Unified  Flow  Solver  with  adaptive  mesh  and 
algorithm  refinement  based  on  direct  numerical  solution  of  the  Boltzmann  equation  coupled  to  kinetic 
schemes  of  gas  dynamics.  Our  strategy  allowed  easy  coupling  of  the  continuum  and  Boltzmann 
solvers  in  a  hybrid  code  with  automatic  domain  decomposition.  We  have  demonstrated  the  UFS 
capabilities  for  several  one-component  gas  flows  and  confirmed  that  our  hybrid  method  results  in 
significant  savings  by  limiting  expensive  kinetic  solutions  only  to  the  regions  where  they  are  needed. 
We  have  shown  that  UFS  can  automatically  introduce  or  remove  kinetic  patches  to  maximize 
accuracy  and  efficiency  of  simulations.  We  have  extended  the  kinetic  solver  to  molecular  gases  with 
rotational  and  vibrational  degrees  of  freedom.  We  have  also  extended  the  UFS  components  to  reactive 
gas  mixtures.  With  some  extra  efforts,  it  seems  feasible  to  produce  an  efficient  solver  for  unified 
simulation  of  practical  problems  of  polyatomic  gas  mixtures  of  different  degrees  of  rarefaction. 
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